深度学习在蛋白质结构预测中的应用及启示

位置特异性打分矩阵(position-specific scoring matrix, PSSM)或位置权重矩阵(position weight matrix, PWM)^[17]是蛋白质及生物信息学里非常重要的统计量. 它主要衡量了不同氨基酸(或核酸)在蛋白质(或DNA)上某个特定序列位置上出现的概率. 在一些机器学习预测蛋白质的二级结构类型时^[18~25]，常会将PSSM作为网络的输入. 但注意PSSM只包含残基绝对位置属性的信息，不包含不同残基配对关联信息.

图1(a)以DNA为例，给出了统计PSSM矩阵的示意流程. 首先给定一个序列库(例如针对基因库的所有DNA数据，或蛋白质库里所有可能的序列)，图中给出了由10个假想DNA序列组成的DNA库；然后统计不同的核酸在特定位置出现的频次矩阵(position frequency matrix, PFM)；再根据PFM得到位置概率矩阵(position probability matrix, PPM)；最后根据图中公式算出位置权重矩阵PWM.

  

Fig. 1  Illustration of (a) position-specific scoring matrix (PSSM) and (b) multiple sequence alignment (MSA). (a) In this illustrative example, PSSM is computed using the formula given in the top right corner based on a DNA database consisting of ten DNA sequences. (b) Sequence alignment (SA) is trying to match the fragment pairs from the two given sequences as much as possible. In the alignment, inserting gaps "-" is allowed. Multiple sequence alignment (MSA) is SA on multiple sequences.

目前大多蛋白质结构预测的深度学习算法的输入中都有多重序列比对信息(multiple sequence alignment, MSA)^{[12,15,26~39]}.

序列比对(sequence alignment)主要任务是针对查询序列(query sequence)从数据库中，用基因信息学的方法找到进化树上尽可能同源的序列，然后根据变异的氨基酸的相似程度，按照特定规则来给该序列与查询序列的相似度打分.

某个序列的变异包括对序列中特定片段的插入、删除和替换. 相对于查询序列，当库里的蛋白质序列变异很少时，则两者相似度高.

当变异多时，还需根据进化同源的特点分类对变异片段进行进一步分析. 变异的氨基酸片段可分为保守片段(功能及化学特性相同)、半保守片段(功能及化学特性相近)和非保守片段(化学特性相差甚远). 显然，若保守片段越多，表明与查询序列越接近.

比对的目标是通过恰当地插入空片段(gap)，使得插入空片段后的2个序列尽量相似(如图1(b)左图所示). 比对的方法有许多^[27]，例如动态规划(dynamic programming)和点阵法(dot-matrix method).

用上述比对方法对若干个给定的序列与查询序列进行比对就称为多重序列比对(multiple sequence alignment, MSA). 通常可用软件ClustalW, MAFFT, ClustalOmega以及MUSCLE等算法程序对多个序列进行MSA比对^[40~52].

而在蛋白质预测中，通常会针对输入的蛋白质序列，从蛋白质数据库中找到与给定序列相近的若干个序列，然后再将这些MSA作为神经网络的输入. 此信息相比于PSSM包含了更为丰富的信息. 可从MSA中看出目标序列大致从哪些序列变异而来. 在深度学习中，MSA数据维度为(N_seq, N_res, 21)，其中N_seq为MSA包含序列的数目，N_res为目标序列的长度，21用于分辨20种氨基酸和gap“-”的热点表征(有时可能为22或23).

如图2所示，图2(b)与2(c)是一个HP蛋白质模型结构^[39]的接触图(Contact Map)与距离图(Distogram). 其中接触图中只有2个残基接触时，才有值(黑)；而Distogram灰度值对应于两残基的距离，当距离大于截断阈值时，灰度为0(白色).

  

Fig. 2  Illustration of contact map and distogram. (a) A typical structure of a given HP protein. (b) The contact of the (c) structure where the black square indicates the matrix element corresponding to two contact residues. (c) The distogram of the (a) structure where the greyness indicates the distance between two residues.

同一序列中不同残基间的接触与否或距离是非常重要的信息，它基本蕴含了蛋白质骨架的三维结构所有的信息. 而且这个信息相比于纯粹的结构三维坐标信息有2个优势：(1) 具有旋转平移不变性，而三维坐标会随着蛋白质的旋转或平移而改变；(2) 表达更简洁及更易标准化. 因为存在关联变异(correlated mutation)现象，有些接触的两氨基酸会同时变异以保证变异后仍接触，故接触图或距离图信息就显得相当重要^[35,36,39].

基于上述原因，在最近的深度学习预测蛋白质结构的实践中^{[12,15,53~62]}，大多都会采用此信息去提高预测准确性或预测给定蛋白质的Contact Map或Distogram.

目前最著名的蛋白质数据库为PDB^[63],即Protein Data Bank，收藏了约1×10⁵多条蛋白质的三维结构数据. 这些结构由X射线、NMR或电子显微镜等方法获得.

Critical Assessment of Protein Structure Prediction (CASP)^[64,65]是蛋白质结构预测科学共同体举办的两年一次的竞赛，每次竞赛优胜者的水平基本代表了当前世界结构预测的最高水准(benchmark progress). 在每次竞赛中，举办方会给出若干个已知结构但未曾公开的蛋白质序列，参赛团队在规定时间内提交各自的结构预测结果，同时不限制预测方法. CASP自1994年以来共举办了14届，其中最近2届的第一名皆来自deepmind的AlphaFold算法.

之前，人们通常用距离均方差root mean squared deviation (RMSD)衡量2个分子构象的接近程度. 但现在模版建模得分template modelling score被认为是更准确的衡量方式^[66]. 其表达式如下：

其中

式中n为蛋白质的残基数，M为旋转平移矩阵. 上式表达的含义是将预测得到的结构与各种旋转平移操作后的真实结构进行比较，取最相近(极大)的那个作为最后的分值.

显然TM score在0~1之间，分数越高表明越准确. 通常认为当TM>0.5时，预测与真实之间的折叠基本一致^[56]；而对同一蛋白质，NMR与X射线测出结构之间的TM分数为0.807 ± 0.107左右. 所以，可认为当TM分数>0.8时，预测的结果已经完全正确.

而AlphaFold2(AF2)近2/3的预测结果达到中低分辨率的实验精度^[12]. 也即AF2几乎解决了单域蛋白质折叠预测问题^[14].

由于多域蛋白质各功能域之间可以相对独立地移动旋转，在评估多域蛋白质结构相似性上，局域距离差异性测试(local distance difference test)是一个比TM分数更佳的评分方式. lDDT不同于TM，不依赖于骨架α碳原子的重叠，能够不受功能域间位移的影响，更加有效地评估结构之间的局域相似性^[67].

本小节仅罗列结构预测涉及到的深度学习技术及原理, 具体请参考相关文献^[68].

神经网络(neural network)可抽象成一个函数y=f(x;w)，它关联了2组信息数据x与y(比如蛋白质的序列x及其结构y)，分别称为网络输入与输出；w为网络的参数. 神经网络训练的目标是为了找到恰当的w使得网络能够根据x准确地预测y.

普适近似原理(universal approximation theorem)^[69]表明单隐藏层的神经网络，只要其激活函数为非线性且神经元数目足够多，便可无限精确近似任意非线性映射. 普适近似原理表明NN可用于拟合任意未知关联.

神经网络设计要点：考察待预测的量y与哪些量有关联，即找出哪些信息可足够推导出y，然后将这些信息与y之间架接合适的神经网络便可. 信息间的关联如果能用现有知识进行关联就用现有知识将其关联；未知关联用神经网络代替.

神经网络选择需要考虑输入输出信息数据特点，目前结构预测中常用的网络结构主要有下面几种.

残差网络(resnet)^[70]的基本思想是不断地将未处理过的信息直接复制并叠加到下面几层由网络抽取出的特征上去. 残差网络于2015年提出，后来被广泛运用于图像处理中.

基于自注意力机制的transformer^[28]近几年备受人工智能领域喜爱，它几乎完全取代循环神经网络^[68]，其基本思想是从不同位置对之间提取信息，适合处理文本类、时序性的信息，不过近年也常用于图像处理. AlphaFold2^[12]中大量使用了自注意力机制.

TBM方法^[14,29~38]大致步骤如图3所示，通常可分为以下几步. 第一步，通过数据库检索，得到目标蛋白质的一组同源性序列(MSA)，并根据MSA获得1个或多个折叠结构模板. 第二步，比对目标序列和模板对应序列，两序列一致的片段直接使用模板的对应折叠结构. 第三步，对于目标序列与模板对应序列不一致的区域，采用碎片组装、优化算法或是数据库方法等单独预测. 当然通常最后还会用诸如分子动力学的优化方法进行模型的精细化(model refinement)，以优化全局结构^[13,14]. 历史上，TBM方法^[14]可以细分成comparative modelling (CM)和threading 2种方法^[29~38]. 其中在CM中，模板与目标序列的同源性较近.

  

Fig. 3  Illustration of basic strategies of template-based method (TBM) and template-free method (FM).

无模板方法^[72~77]的流程见图3右半侧，FM从蛋白质数据库中依MSA比对结果找到一些片段的结构并将其放入片段库中，然后找到评分较高的片段结构拼成初始结构^[72~74]，接着采用FM里非常重要的片段组装(fragment assembly)方法^[72]，大致冻结片段的结构并以片段结构为单元来演化全局结构，比如可根据粗粒化的势能函数用梯度下降化进行能量优化.

人们发现在蛋白质变异过程中经常出现关联变异(correlated mutation)的现象：一条蛋白质链内若发生变异，总是2个氨基酸成对地变异；因为演化压力会迫使蛋白质维持一致构型，原本接触的氨基酸对在变异过程中继续保持接触，可以避免其形状发生剧烈变化. 因此，这就使得残基接触对(inter-residue contact map)的信息极为重要^[56~61].

早期有许多传统方法致力于预测残基接触对. 处理该问题的早期算法，倾向于以一次一对的形式、孤立地预测每个接触对是否可能. 由于忽视了蛋白质包含的全局信息：一个残基对是否接触受到序列中其他残基的影响，早期算法陷入了困境，预测效果糟糕. 而之后研究者提出了充分利用全局信息的预测方法，例如基于Markov随机场模型MRF的direct coupling method (DCA)^[58~61]，在残基接触预测上获得了突破性的成就.

深度神经网络在预测残基接触对问题上，也表现出了异常优异的性能，有时甚至还直接被用于预测键角等信息. 这些预测特征均可作为约束，辅助指导无模板方法.

比如，Raptor X-Contact深度学习模型^[39]将Contact Map的预测当成图片分割任务来对待，Raptor X-Contact所采用的方法也被其他方法，如ResPRE^[54]所采纳. ResPRE^[54]采用了图片识别领域非常著名的残差网络(Resnet)模块^[70]，残差网络的重要思想是不断地将网络前面的信息直接复制到网络后面.

而AlphaFold1^[15]又将Contact Map拓展成距离直方图(distogram)预测，基于此，它在2018年CASP13的比赛中获得了巨大成功.

2020年的CASP14的比赛中，AlphaFold2 (AF2)^[12]取得了骄人的成绩. 对来自89个域(domain)实验测得的蛋白质结构，AlphFold2在88个域TM分数>0.5，59个域分数>0.914. 前者意味着预测结果与答案之间折叠基本一致. NMR、X射线晶体学测出的一组112个单域蛋白质，序列相同率大于95%. NMR与X射线测出的结构之间的TM值为0.807±0.107. 这说明AlphFold2的近60%的预测达到中低分辨率的实验精度. 也就是说AlphFold2几乎解决了单域蛋白质折叠预测问题^[14].

AlphaFold2深度学习模型的结构简图如图4所示，具体参考文献^[12]. 它分别借助了基因同源信息和蛋白质结构数据库模板信息. 如图所示，根据同源信息，可得到序列比对信息MSA，通过同源搜索得到与输入序列同源相近的(s-1)条序列和输入序列一起放到MSA数组里，再通过线性神经网络变换得到MSA表征，此表征的维度为(s，r，c)，其中r为蛋白质的序列长度，c为表征的特征数(通道数). MSA表征包含了输入序列与其他同源序列间的关系.

  

Fig. 4  Sketch of the AlphaFold2 model. Detailed description is referred to Ref.[12].

而另一输入通道中，主要输入与MSA相对应的序列的结构残基对距离信息以及扭转角的信息. 在具体输入时，AF2将距离对长度划分成64个离散块(64 bins)，并将其转化为概率的形式，故对应数组形状为(s，r，r，64)，取值为0~1. 注意配对表征中，只包含了MSA除输入序列之外的某个序列自己结构信息，不同序列之间并没有进行信息的关联.

然后再将MSA表征与配对表征输入一个称为Evoformer的模块，此模块主要将MSA的信息(同源性差异)与结构信息整合起来，最后得到输入序列的MSA表征与输入序列的配对表征. 此时，输入序列的配对表征同时将演化信息与其他模板结构信息有机地融合在了一起. Evoformer主要利用了自注意力机制来实现上述信息整合.

而下一个结构模块structure module主要的功能是将Evoformer预测的配对表征展开成三维空间结构，同时亦承担一定的预测调整功能. 此模块的结构大致如图5所示. 一条蛋白质骨架结构可想象成一系列三角形的叠加，三角形的中心相当各个残基α碳的坐标，三角形平面本身代表N-α-C-C构成的三角形. 这样，此骨架可由2个数组表示，数组形状分别为(r，3×3)和(r，3)，分别表示每个三角形取向与位置.

  

Fig. 5  Illustration of how the pairing information is transformed into the 3D structure using neural networks in AlphaFold2^[12].

初始时，假设所有氨基酸都在原点，然后将此初始骨架与配对表征输入结构模块，由于配对表征存有距离对及取向信息，故可通过一个称为不变点注意力神经网络模块将其初步还原成展开的骨架结构，紧接着再加入侧链原子从而得到全原子的三维结构.

如图4所示，最后再将中间输出的MSA信息、配对信息和3D结构信息重新叠加输入到Evoformer，如此反复迭代3次，最终到预测结果.

因为PDB中只有大约1×10⁵多个的序列有对应的三维结构数据. 而在big fantastic database (BFD)蛋白质序列数据有多达2,204,359,010个序列，虽然这些序列并不一定有对应的三维结构信息(无标签)，但self-distillation dataset的训练技巧可以将这些无答案的题目作为作业进行训练，自己提高预测准确度，AlphaFold2用此扩大训练集并进一步提高了预测准确度.

后来有诸多研究团队对AlphaFold2进行了拓展与提升. 例如：Baker团队^[78]的RoseTTAFold发展了三通路神经网络(three-track neural network)，对AlphaFold2只包括1D序列信息和2D距离图信息的两通路神经网络模型进行了拓展，引入了3D结构通路道网络模块；高毅勤团队的MindSpore算法^[79]对AlphaFold2的计算速度进行了较大的提升.

真实蛋白质结构预测无论从训练数据准备还是模型构建及训练都极其复杂. 因此，人们希望找一个简单的蛋白质模型，以便能快速地试验他们的想法. 就如手写数字识别(对应数据集为MNIST)^[80]对于图像识别一样，所有的方法都会用MNIST数据集先来检验其有效性.

而HP蛋白质模型就是这样的模型^[2,6,8,9]. 它仅有2类氨基酸H和P，其中H代表厌水型氨基酸，P代表亲水型.

我们基于此HP模型，提出了一个强关联神经网络^[16]，如图6所示，此神经网络有2个核心要素，一是不同于传统的向量表征，它采用一个小的神经网络来代表每个氨基酸，每种氨基酸都用一个神经网络来表征，不同氨基酸对应的网络的权重亦不同，而相同的氨基酸共享网络权重；二是它有一个自洽循环通路，这样可使得输出的信息(环境)与氨基酸的属性发生强关联.

  

Fig. 6  Architecture of the strongly-correlated neural network (SCN) where r indicates number of residues and c indicates number of features or channels.

该研究发现与传统向量表征方法相比，强关联网络极大提升了预测准确性，提高了约20个百分点.