AlphaFold2怎么给蛋白结构打分？
前言：
六大派围攻光明顶，囧蒋迫孤身战群英。
上个月的通稿一度密集得让人应接不暇。
12.8 深势科技发布 Uni-Fold，在CASP14的目标集上，平均 Ca-LDDT 打分 82.6；
12.9 天壤智能发布 TRFold2，在CASP14的目标集上，平均 TM-score 打分 82.7；
12.14 华深智药发布 HeliXonAI，在最新的CAMEO测试集上，LDDT和TM-score两项打分均超过 AlphaFold2。
除自媒体，纸媒也有大篇幅报道。《人民日报》在 12.27 第18版报道了天壤智能《人工智能预测蛋白质结构》，将通稿推向了新的高度。
这里必须澄清：
LDDT 的分值介于 0 ~ 100，分值越高越好。
TM-score 的分值介乎于 0 ~ 1，也是分值越高越好，将分值 x 100 也行吧。
如此一来，深势 Uni-Fold 和天壤 TRFold2 不能直接比较高下，两种打分无法直接比较。不能说因为天壤 TRFold2 的82.7高于深势 Uni-Fold 的82.6，所以前者优于后者。
考虑到天壤的通稿晚深势一天发布，笔者不能不认为这样巧合的打分才是其将 TM-score x 100 的原因。
笔者仔细查了文献中报道的 AlphaFold2 的得分。当然，不必查文献。CASP 官网保留有所有参赛队伍的每一个模型的预测结果，不怕麻烦，自己下载文件算一下也行。
AlphaFold2 的得分如下：
LDDT-Ca 的中位数为 90.3，平均值 82.6；[1] 补充材料 P53 RMSD95 的中位数为 0.96 Å（主链），平均值 1.5 Å（全原子）；[1] 正文 P2 全链 TM-score 平均值为 0.91；[2] 正文 P7 结构域 GDT_TS 中位数为 92.4，平均值 87.66；[3] 摘要，[1] 补充材料 P46 加和 z-score (>2.0) 为 244.0；[3] 摘要
由此可见:
使用 LDDT-Ca 作为结构测度，得分同为 82.6，深势 Uni-Fold 已经媲美 AlphaFold2；考虑其运算速度是 AlphaFold2 的 2 ~ 3倍，可以认为 Uni-Fold 相较 AlphaFold2 已略胜一筹。并且，深势还将 Uni-Fold 的训练代码开源了 —— 目前市面上唯一开源的基于 AlphaFold2 的训练模型，大赞。
使用 TM-score 作为结构测度，天壤 TRFold2 得分 0.827，仍然低于 AlphaFold2 的 0.91，略高于 RoseTTAFold 的0.813，因此可以认为 TRFold2 仍稍逊 AlphaFold2。
华深 HeliXonAI 则直接给出 LDDT 和 TM-score 两项得分，并测试了开源的深势 Uni-Fold；在 CAMEO 官方给出的评分下，其以不小的幅度超越 AlphaFold2，可谓后来居上。
笔者想澄清的第二点是，千万不要以为“创新改进”不容易，而“复现”就很简单。算法设计、代码实现、程序部署，每一步都不简单；何况“复现”需要真金白银、积月累日的投入。另一方面，单单“复现”也有培养训练人材的作用，“唯手熟尔”就很有意义。
第三，依据各家公司的新闻稿，综合排名当下世界的蛋白质结构算法：
HeliXonAI > Uni-Fold ≥ AlphaFold2 >
TRFold2 ≥ RoseTTAFold > 其它各种 fold
当然，除了华深、深势、天壤公布了各自的蛋白结构预测的训练模型，国内还有多支队伍正在开发自己的算法。笔者期待这些队伍在今年的 CASP15 中大放异彩。
以上为引子。
笔者想借此机会说一下蛋白质结构的打分函数，或者说，蛋白质结构相似度 / 差异度的度量。
严格意义上，本文不是对 AlphaFold2 的拆解，因为这些打分函数本已经广为应用，不为 AlphaFold2 独创独有。但讲解这些函数有意义，因为 AlphaFold 团队的确以此为基础做出了独特而有效的创新，比如 pLDDT。
并且，从历史革沿看，笔者想带大家读读老文献，让学习者更能理解这些函数之所以被保留至今的原因。They’re kept using for a reason.
正文目录
AlphaFold2 所用的蛋白质结构打分函数 Z-score RMSD95 RMSD RMSD 的改进：rRMSD、Z-rRMSD TM-score LG-score MaxSub TM-score = 0.5 的意义 GDT_TS LGS = LCS + GDT LDDT RMSD 的变体：dRMSD DLDDT pLDDT 信息熵 总结
本节参考文献
[1] Jumper et al., Highly accurate protein structure prediction with AlphaFold, Nature (2021)
[2] Kryshtafovych et al., Critical assessment of methods of protein structure prediction (CASP)—Round XIII, Proteins (2021)
[3] Jumper et al., Applying and improving AlphaFold at CASP14, Proteins (2021)
AlphaFold2 所用的蛋白质结构打分函数
Z-score，也称 standard score，是观测值与平均值之差除以标准差（公式0）。统计学中，Z-score 表征一个观测值多大程度上偏离平均值。
CASP14 用 Z-score 表示各个参赛方法的总结果，这其实并非评估预测结构与天然结构相似几何。—— AlphaFold2 以 244.0 的 Z-score 加和总分遥遥领先。这意味着，其几乎每个预测结果都远远好于所有参赛队伍的平均结果。
RMSD 是最常用、最经久不衰的蛋白结构度量，也最被诟病。它定义简单，就是两个结构对应原子坐标差值的均方根（公式1）。
定义分子的“惯性半径”（Rg，radius of gyration）如公式2，那么，RMSD 的公式可以由 Rg 表达（公式3）。这里注意，两个结构以最佳方式重叠后，质心应当重叠，因此可以不失一般性地令其质心在原点，在公式2中约去向量 rc。[4]
因为计算 RMSD 必须先将相比较的2个结构最大化地重叠，即找到一个刚体变换，令 RMSD 最小，所以 RMSD 的计算公式可以更数学化地写为公式4（考虑原子质量 mi）。
AlphaFold2 用的是 RMSD95，即只考虑95%的残基最大化的重叠（公式5），这样降低了一些灵活区域的影响。[1]
虽然 RMSD 计算简单，物理意义明了，但它有一些明显的缺点：
计算强烈依赖于两个结构的重叠对齐算法；
无法区分整体的结构拓扑；
对结构的缺失部分不敏感；
有链长依赖。
对2、3，通常，边角灵活区域的结构浮动所造成的 RMSD 差异，可能盖过整体结构的重叠。也就是说，即便2个结构整体相似（有相同的fold），其二者也会因为少部分灵活区域的不同而有较大的 RMSD。
对4，同样的 RMSD 值对较小的蛋白可能意味着巨大的结构差异，对较大的蛋白则可能代表微小的结构差异。这一点我们可以由公式4看出—— Rg 显然依赖于蛋白的大小。实际上，<Rg> ~ N^0.39。[5]
为了回应 RMSD 的上述缺点，诸多 RMSD 变体被设计出来，我们介绍2个：rRMSD和Z-rRMSD。
2001年，Skolnick 与学生 Betancourt 设计出 relative-RMSD，即 rRMSD [4]，来消除链长依赖 —— 着眼点就是公式3。
考虑“对齐相关系数”（aligned correlation coefficient, ACC），定义一组随机结构的 ACC 的平均值为 C(N)（公式6）。在考察了1300个同源性低于30%的结构的连续残基片段后，作者获得了 C(N) 的一个估值（公式7）。—— 21年后的今天看，这个公式是这篇文章最有价值的成果，而不是它的 rRMSD。
公式7的两个特征幂指数 4.7 和 37 的单位是“残基”。
4.7 个残基大约对应 2.8 nm，与球蛋白的 Lp 值接近（persistence length）。Lp 是高聚分子的分子链灵活性的特征值。它表示一条高分子链上的两点，在相距一定距离后，其各自的取向彼此不相关（decorrelated）。设想一手拽住一根长绳抖动，绳上相距很近的位置的取向必然相互关联，换句话说，某一处的位置和运动速度，将在一定距离内，影响它后面部分的位置和运动；但是，当两点相距足够远，这种相关性就没有了。
37 个残基，作者认为，是其所考察的中等长度多肽链（片段）的相关长度。虽然有强行解释的味道，但从物理含义的角度理解公式依然很有意思。
有了C(N)，可以估计随机的两个多肽结构可能有的平均 RMSD（公式8），进而定义 rRMSD（公式9），这样就消除了链长依赖，改进了 RMSD 的缺点4。
一脉相承地，2004年，Skolnick 与学生张阳基于 rRMSD，设计了 Z-rRMSD（公式10，公式11），一种类 Z-score 的统计重要性（statistical significance）的度量。[5]
Z-rRMSD 从概率上判断2个结构（片段）是否相似，更重要地，是否基于已知所有结构”universally”相似。如果2个结构的 Z-rRMSD 打分很高，那么这2个结构很可能有相同的 fold。也就是说，Z-rRMSD 改进了RMSD 的缺点2、3。
”有相同的 fold“？这不是 TM-score 所擅长的吗？对。文献5就是提出 TM-score 的著名论文。Z-rRMSD 可以作为 TM-score 的替代。
本节参考文献
[4] Betancourt & Skolnick, Universal similarity measure for comparing protein structures, Biopolymers (2001)
[5] Zhang & Skolnick, Scoring function for automated assessment of protein structure template quality, Proteins (2004)
很少有研究是完全新颖的，TM-score 也不例外，它其实是LG-score的改进MaxSub的改进。
1997年，Michael Levitt，这位日后的诺贝尔奖得主，与他的博后Mark Gerstein，目前耶鲁大学生物医学信息学、分子生物物理与生物化学、计算机学的三聘教授，出于统一比较蛋白序列相似度和结构相似度的想法，提出了Levitt-Gerstein score，也就是 LG-score（公式12）。[6] LG-score 是首个不依赖于序列的结构对齐打分。
2000年，Siew等人提出 MaxSub 算法，包括一种结构重叠对齐算法和一种基于 LG-score 的结构相似度打分（公式13）。[7]
2004年，张阳与导师 Skolnick设计出template modeling score，也就是 TM-score（公式14），评估全长模型的预测，消除蛋白质的大小对结构打分的影响；更重要地，TM-score 可判断2个相比较的结构是否属于同一 fold，即在整体结构或拓扑层级评价结构。[5]
相对于 LG-score，MaxSub 去掉了未匹配残基的空位罚分（gap penalty），并且去掉了无用的常数 M，代之以 1/LN，从而在优化过程中需要最大化地重合两个比较结构的子结构，只考虑对齐后 di ≤ d0 的残基，尽可能地增加 LT；相对于 MaxSub，TM-score 则将常数 d0 替换为一个 LN 的函数，从而使得整体打分没有链长依赖。
理解 TM-score 的物理含义以及它相对 LG-score 和 MaxSub 的改进，关键在于理解公式14中的 d0 —— d0 是天然结构的长度（LN）的函数，可是为什么令 di 除以 d0 就能消除链长依赖？
注意公式14中 di 的定义，它其实就是单个残基的 RMSD，它的均值我们在公式8中已经给出。因此，di 的均值可以由公式15近似表示。作者发现，令常数 h = 0.75，则 di 的均值可以进一步用更简单的函数形式近似，即公式14中的 d0。
现在我们清楚了 TM-score 的物理意义：
LG-score + rRMSD —> TM-score
也就是说，依然如 RMSD 一般对齐两个相比较的结构，只不过这两个结构未必含有相同个数的残基，也就是可以只对齐结构的一个子集；对齐后，依然考虑两两对齐的残基之间的距离，只不过将这个距离除以它的统计平均值，并且这个平均值是链长的函数。
那么，由此可以凭直觉猜测 TM-score 有3大优势（改进）：
不依赖链长，具体原因我们上面分析了。
不依赖序列（sequence independent），不需要两个相比较的结构有同样数目的残基。
最大程度地对齐两个结构的子集 —— 可衡量两个结构的整体拓扑。因为必须打分必须要除以 LN，所以，需要尽可能找到两个结构的最大的相似子集。
6年后，2010年，针对上述的第3点，张阳与学生 XU Jinrui 考察了 TM-score 在判断蛋白质 fold 方面的意义，量化了 TM-score 打分的统计意义：TM-score = 0.5 意味着两个比较的结构“很可能”同属一个 fold。[8] 这篇文章可能是 TM-score 获得巨大关注并被广为应用的原因 —— 我把我自己阐释清楚，从而别人也明白我。己知，而后人知。
本节参考文献
[5] Zhang & Skolnick, Scoring function for automated assessment of protein structure template quality, Proteins (2004)
[6] Levitt & Gerstein, A unified statistical framework for sequence comparison and structure comparison, PNAS (1998)
[7] Siew et al., MaxSub: an automated measure for the assessment of protein structure prediction quality (2000)
[8] Xu & Zhang, How significant is a protein structure similarity with TM-score = 0.5? Bioinformatics (2010)
CASP5 之后，Zemla发文详细介绍了2种互为补充的算法：LCS & GDT，统称LGA（local / global alignment）[9]，其中 GDT 早在 CASP3 就被作为 RMSD 的一个改进，目标是开发一种对蛋白的结构拓扑敏感的度量 [10]。
LCS 和 GDT 的定义很简单：
LCS: longest continuous segments under specified RMSD cutoff
GDT: global distance test under specified distance cutoff
LCS 是在预设的 CA-RMSD 阈值下，对齐两个分子的由连续残基组成的子结构，计算对齐的子结构之间的 RMSD，令符合此条件的连续残基片段最长。
可以设想：当这个预设的 RMSD 阈值很大时，两个分子对齐后，二者的 RMSD 小于阈值，那么 LCS 就是整条链；当这个阈值很小时，再如何调整可能也只能使得连续的几个残基在对齐后的 RMSD 小于阈值。极端情况就是2个残基，当然一条线无论如何可以对齐。
GDT 是在预设的 CA 距离阈值下，对齐两个分子，考虑对齐后有多少对残基之间的距离小于阈值，这些残基可以不连续，将符合条件的残基数除以分子的残基总数即可。
类似地，当距离阈值很大时，两个分子对齐后，所有对应的残基两两之间的距离都小于阈值，那么 GDT = 1；当距离阈值很小时，无论如何调整（整个分子的）对齐方式，最终散落分布在整个分子上的、符合距离条件的残基也不会很多，GDT —> 0。
从定义可以看出：LCS 关注的是局部（子）结构；GDT 比较的是全局结构。GDT 可以评估整体结构的拓扑，则难以反映局部结构的差异，如蛋白骨架构象、侧链的堆积等。
计算公式非常简单，其中 GDT_x 表示阈值为 x 的 GDT（公式16）：
但是，想必大家也意识到了，获取 LCS 或者计算 GDT 的关键是“对齐算法”。当然可以用 RMSD 或者 TM-align 所用的对齐算法，再进行多轮迭代。Zemla 也是这么做的。考虑到这是一个 NP-hard 问题，基于 heuristic 的算法不能保证找到最优解。
2011年，李明教授和学生许锦波、卜东波、李帅成共同开发了optGDT算法，证明可以在多项式时间内近乎解析地找到最优解。[11]
本节参考文献
[9] Zemla, LGA: a method for finding 3D similarities in protein structures, Nucl. Acids Res. (2003)
[10] Zemla et al., Processing and analysis of CASP3 protein structure predictions, Proteins (1999)
[11] Li et al., Finding nearly optimal GDT scores, J. Comput. Biol. (2011)
LDDT 的全称是 local distance difference test，局部距离差测试，在它被提出来的时候也不是全新概念，它来自 dRMSD + GDT。
dRMSD（distance RMSD）是 RMSD 的变体，计算每一个分子内部原子对的距离，再考虑2个分子的对应的原子对的距离差 [12]。也就是说，RMSD 计算的是两个分子对应原子的距离差（在重叠对齐后），dRMSD 计算的是两个分子对应原子对的距离差的差值（公式17）。dRMSD 的最大的特点是：它的计算不依赖于将两个结构重叠对齐！
2010年，CASP9，Mariani 等人在评估基于模板的预测结构时引入了 DDT（distance difference test）[13-14]。作者将 DDT 视为 dRMSD equivalent of GDT。
dRMSD + GDT —> LDDT
GDT 是两个分子重叠后小于预设距离阈值的残基数与残基总数的比值。DDT 是小于预设距离阈值的分子内距离对的数目与距离对总数的比值。
Mariani 等人在2013年正式推出 LDDT 的论文中并没有给出数学公式，仅有文字描述，甚至不是伪代码算法；Senior 等人在 AlphaFold1 的论文的补充材料给出了一个公式（公式18）[15]。因为 AlphaFold1 是在王晟和导师许锦波的以接触图+ResNet的框架下的改进，而其改进之一是用距离图（distogram）取代接触图，所以仿照 LDDT 定义了 distogram LDDT，即 DLDDT（公式19）。
很快，依旧沿着王-许方法的路径，杨建益和导师 Baker 用取向图（orientogram）取代了距离图；紧接着，CASP14 中，AlphaFold2 另辟蹊径，使用了完全不同的深度学习框架。可以预见，distogram 与 DLDDT 都将与过去的许多度量一样，成为历史名词。
就公式18我们解释一下 LDDT 的含义：
考虑天然结构内的原子之间的距离 Dij，需要 Dij > 15 Å 且 |i - j| ≥ r，由此得到符合条件的原子对 (i, j)。这时考虑预测结构内这些原子对的距离 dij，如果 |Dij - dij| < tolerance，则加1。分别令 tolerance = 0.5, 1, 2, 4 Å，计算满足条件的距离差的比例，再取均值。
如果 tolerance 很小，那么意味着满足 |Dij - dij| < tolerance 的原子对 (i, j) 在天然结构和预测结构中的距离几乎一致。
从另一个角度看，无论 tolerance 预设为何值，满足条件的原子对 (i, j) 越多，说明在此容忍度下，预测结构与天然结构内部原子之间的距离越一致，那么两个结构应当越相似。—— 通过计算预测结构的内禀距离判别是否与天然结构相似。
由此，在 AlphaFold2 中，Jumper 等人设计了pLDDT，per-residue LDDT，通过计算预测结构的内禀距离，在缺乏天然结构的情况下，判断预测结构中的每一个残基的可信度。[1]
利用 predicted pLDDT 进行预测精度的自评估是 AphaFold2 正文列举的7大创新之一。
实际上，细读 AlphaFold2 论文可以看到，predicted pLDDT 并非作者的首选。—— 信息熵才是。信息熵的基本思想就是比较两个结构的内禀距离的概率分布的差异，与 LDDT 的基本想法是一样的。
本节参考文献
[12] Bordogna et al., Predicting the accuracy of protein-ligand docking on homology models, J. Comput. Chem. (2010)
[13] Mariani et al., Assessment of template based protein structure predictions in CASP9, Proteins (2011)
[14] Mariani et al., lDDT: a local superposition-free score for comparing protein structures and models using distance difference tests, Bioinform. (2013)
[15] Senior et al., Improved protein structure prediction using potentials from deep learning, Nature (2020)
总结
总结 AlphaFold2 所用的结构度量：
RMSD95，CASP1 开始使用，比较整体结构，需要重叠对齐。 GDT_TS，CASP3 开始使用，比较整体结构，需要重叠对齐，判断 domain 相似度。 TM-score，CASP5 开始使用，比较整体结构，需要重叠对齐，可判断整条链是否同属一个 fold。 LDDT，CASP9 开始使用，比较局部结构，不需重叠对齐，预测值可自我评估预测可信度。
我们特意提到每一种度量被引入使用的时间。例子还可以多举一些（当然是不完全列举）：
CASP1 (1994): RMSD CASP3 (1998): GDT MaxSub CASP5 (2002): TM-score CASP7 (2006): HBscore - 1st score based on atomic interactions CASP8 (2008): MolProbity MCRS: main chain reality score HBmc & HBsc: hydrogen bond correctness corRot: rotates correctness GDC-sc: side chain positioning CASP9 (2010): LDDT CAD score SphereGrinder score CASP10 (2012): RPF / DP-score
正如我们在上面提到的注定今后不再有用的 CASP13 中的 DLDDT，上面这许多打分，一方面应预测算法的精度提高而被提出，用于考量预测结构的不同方面（fold, domain, local, etc.），另一方面随算法的迭代和进步而被快速地抛弃。
事实上，上面列举的一些结构度量毋宁说是“反应坐标”。例如，HBscore 在论文（通稿）中被宣布为有史以来第一个能反映原子间相互作用的结构打分。其实就看天然结构中存在的氢键对是否也存在于预测结构中，二者数目的比值就是 HBscore。笔者甚至在一篇文章的引言论述 GDT 等打分的缺陷时说其不能反映原子间相互作用。怪无语的。
一直保留至今的结构打分函数，笔者认为满足了三项指标，简称CSP 原则：
计算经济（Computationally efficient），
区分度高（Structurally differentiable），
物理意义明显（Physically meaningful）。
大部分结构度量的计算都很快；少数，如信息熵，因为需要计算一定区域内所有原子对之间的距离，因而时间复杂度 ≥ O(N^2)。这一类度量对较大体系，特别当需要将此度量应用于动力学轨迹时，计算较慢。蛋白质的结构打分的计算最费时的部分通常都体现在重叠对齐算法上；因为其是 NP-hard，所以一般用 heuristic 解决。
区分度高毋庸多言，它有两方面含义：（1）其打分与多数主流打分函数呈线性关系，也就是说，不至于在 RMSD 看来一个很糟糕的结构和一个很好的结构，用这种打分的结果竟然非常接近；（2）能够或者在局域区分细小的结构变化，或者在整体区分拓扑变化。
最后，物理意义明显，也就是符合物理直觉。显然，定义越复杂的度量，越远离直觉，尽管可能越有用。
本节参考文献
[16] Murzin et al., SCOP: a structural classification of protein database for the investigation of sequences and structures
末：善战者无赫赫之功。
本文完。
2022.1.9 于深圳
AlphaFold
AlphaFold2 硬核技术
蛋白质科学史