天然蛋白质通过调节一维氨基酸序列信息,能够精准地制备具有特殊的三维空间结构的蛋白质分子,实现特定的生理功能.而蛋白质结构预测希望代替大自然通过各种方法从一维序列信息推断其三维空间结构.蛋白质结构预测问题提出至今已困扰我们五十多年[1~3].
通常认为蛋白质折叠的驱动力包括以下几种[2~5]:氢键作用、分子间的范德华相互作用、残基骨架扭转角的选择性、静电作用、非极性基团的厌水相互作用和构象熵.上述驱动力可被统一地描述为“力场”或势能函数.此势能函数也被称为蛋白质折叠能量全景图(protein-foldingenergylandscape).而统计热力学研究表明,此全景图呈漏斗形[2,6~10].大部分未折叠构象形成了高能量地势较缓的平原;而少数折叠构象形成能量低且地势陡峭的漏斗底部.
理解了快速折叠的原理不代表解决了蛋白质预测问题.
在传统的蛋白质折叠预测中,人们通常经过构造或选择力场,从某非天然态出发,用各种动力学计算或模拟方法(例如分子动力学模拟)演化其构象,直至能量达到全局最小[1~3].但传统预测方法会随着残基数目增加计算量迅速上升,事实上传统方法对大多蛋白质结构预测都无能为力[12].
此困境一度让蛋白质折叠预测领域的人们绝望.因此,人们不再依赖基于纯粹物理机制的方法,而是采用结合数据驱动的方式[13,14].最近十多年,这种结合数据驱动的方法随着深度学习在2012年的兴起而愈受重视.直至近3年,AlphaFold[12,15]的突然崛起,特别是AlphaFold2预测蛋白质的高准确性甚至让许多人相信蛋白质折叠预测难题将被解决[14].
位置特异性打分矩阵(position-specificscoringmatrix,PSSM)或位置权重矩阵(positionweightmatrix,PWM)[17]是蛋白质及生物信息学里非常重要的统计量.它主要衡量了不同氨基酸(或核酸)在蛋白质(或DNA)上某个特定序列位置上出现的概率.在一些机器学习预测蛋白质的二级结构类型时[18~25],常会将PSSM作为网络的输入.但注意PSSM只包含残基绝对位置属性的信息,不包含不同残基配对关联信息.
图1(a)以DNA为例,给出了统计PSSM矩阵的示意流程.首先给定一个序列库(例如针对基因库的所有DNA数据,或蛋白质库里所有可能的序列),图中给出了由10个假想DNA序列组成的DNA库;然后统计不同的核酸在特定位置出现的频次矩阵(positionfrequencymatrix,PFM);再根据PFM得到位置概率矩阵(positionprobabilitymatrix,PPM);最后根据图中公式算出位置权重矩阵PWM.
Fig.1Illustrationof(a)position-specificscoringmatrix(PSSM)and(b)multiplesequencealignment(MSA).(a)Inthisillustrativeexample,PSSMiscomputedusingtheformulagiveninthetoprightcornerbasedonaDNAdatabaseconsistingoftenDNAsequences.(b)Sequencealignment(SA)istryingtomatchthefragmentpairsfromthetwogivensequencesasmuchaspossible.Inthealignment,insertinggaps"-"isallowed.Multiplesequencealignment(MSA)isSAonmultiplesequences.
目前大多蛋白质结构预测的深度学习算法的输入中都有多重序列比对信息(multiplesequencealignment,MSA)[12,15,26~39].
序列比对(sequencealignment)主要任务是针对查询序列(querysequence)从数据库中,用基因信息学的方法找到进化树上尽可能同源的序列,然后根据变异的氨基酸的相似程度,按照特定规则来给该序列与查询序列的相似度打分.
某个序列的变异包括对序列中特定片段的插入、删除和替换.相对于查询序列,当库里的蛋白质序列变异很少时,则两者相似度高.
当变异多时,还需根据进化同源的特点分类对变异片段进行进一步分析.变异的氨基酸片段可分为保守片段(功能及化学特性相同)、半保守片段(功能及化学特性相近)和非保守片段(化学特性相差甚远).显然,若保守片段越多,表明与查询序列越接近.
比对的目标是通过恰当地插入空片段(gap),使得插入空片段后的2个序列尽量相似(如图1(b)左图所示).比对的方法有许多[27],例如动态规划(dynamicprogramming)和点阵法(dot-matrixmethod).
用上述比对方法对若干个给定的序列与查询序列进行比对就称为多重序列比对(multiplesequencealignment,MSA).通常可用软件ClustalW,MAFFT,ClustalOmega以及MUSCLE等算法程序对多个序列进行MSA比对[40~52].
而在蛋白质预测中,通常会针对输入的蛋白质序列,从蛋白质数据库中找到与给定序列相近的若干个序列,然后再将这些MSA作为神经网络的输入.此信息相比于PSSM包含了更为丰富的信息.可从MSA中看出目标序列大致从哪些序列变异而来.在深度学习中,MSA数据维度为(Nseq,Nres,21),其中Nseq为MSA包含序列的数目,Nres为目标序列的长度,21用于分辨20种氨基酸和gap“-”的热点表征(有时可能为22或23).
如图2所示,图2(b)与2(c)是一个HP蛋白质模型结构[39]的接触图(ContactMap)与距离图(Distogram).其中接触图中只有2个残基接触时,才有值(黑);而Distogram灰度值对应于两残基的距离,当距离大于截断阈值时,灰度为0(白色).
Fig.2Illustrationofcontactmapanddistogram.(a)AtypicalstructureofagivenHPprotein.(b)Thecontactofthe(c)structurewheretheblacksquareindicatesthematrixelementcorrespondingtotwocontactresidues.(c)Thedistogramofthe(a)structurewherethegreynessindicatesthedistancebetweentworesidues.
同一序列中不同残基间的接触与否或距离是非常重要的信息,它基本蕴含了蛋白质骨架的三维结构所有的信息.而且这个信息相比于纯粹的结构三维坐标信息有2个优势:(1)具有旋转平移不变性,而三维坐标会随着蛋白质的旋转或平移而改变;(2)表达更简洁及更易标准化.因为存在关联变异(correlatedmutation)现象,有些接触的两氨基酸会同时变异以保证变异后仍接触,故接触图或距离图信息就显得相当重要[35,36,39].
基于上述原因,在最近的深度学习预测蛋白质结构的实践中[12,15,53~62],大多都会采用此信息去提高预测准确性或预测给定蛋白质的ContactMap或Distogram.
目前最著名的蛋白质数据库为PDB[63],即ProteinDataBank,收藏了约1×105多条蛋白质的三维结构数据.这些结构由X射线、NMR或电子显微镜等方法获得.
之前,人们通常用距离均方差rootmeansquareddeviation(RMSD)衡量2个分子构象的接近程度.但现在模版建模得分templatemodellingscore被认为是更准确的衡量方式[66].其表达式如下:
其中
式中n为蛋白质的残基数,M为旋转平移矩阵.上式表达的含义是将预测得到的结构与各种旋转平移操作后的真实结构进行比较,取最相近(极大)的那个作为最后的分值.
显然TMscore在0~1之间,分数越高表明越准确.通常认为当TM>0.5时,预测与真实之间的折叠基本一致[56];而对同一蛋白质,NMR与X射线测出结构之间的TM分数为0.807±0.107左右.所以,可认为当TM分数>0.8时,预测的结果已经完全正确.
而AlphaFold2(AF2)近2/3的预测结果达到中低分辨率的实验精度[12].也即AF2几乎解决了单域蛋白质折叠预测问题[14].
由于多域蛋白质各功能域之间可以相对独立地移动旋转,在评估多域蛋白质结构相似性上,局域距离差异性测试(localdistancedifferencetest)是一个比TM分数更佳的评分方式.lDDT不同于TM,不依赖于骨架α碳原子的重叠,能够不受功能域间位移的影响,更加有效地评估结构之间的局域相似性[67].
普适近似原理(universalapproximationtheorem)[69]表明单隐藏层的神经网络,只要其激活函数为非线性且神经元数目足够多,便可无限精确近似任意非线性映射.普适近似原理表明NN可用于拟合任意未知关联.
神经网络设计要点:考察待预测的量y与哪些量有关联,即找出哪些信息可足够推导出y,然后将这些信息与y之间架接合适的神经网络便可.信息间的关联如果能用现有知识进行关联就用现有知识将其关联;未知关联用神经网络代替.
神经网络选择需要考虑输入输出信息数据特点,目前结构预测中常用的网络结构主要有下面几种.
残差网络(resnet)[70]的基本思想是不断地将未处理过的信息直接复制并叠加到下面几层由网络抽取出的特征上去.残差网络于2015年提出,后来被广泛运用于图像处理中.
基于自注意力机制的transformer[28]近几年备受人工智能领域喜爱,它几乎完全取代循环神经网络[68],其基本思想是从不同位置对之间提取信息,适合处理文本类、时序性的信息,不过近年也常用于图像处理.AlphaFold2[12]中大量使用了自注意力机制.
传统的蛋白质结构预测方法[14]主要基于以下2种模型:基于模板的方法(template-basedmethod,TBM)[29~38]和无模板方法(template-freemethod,TFM)[72~77];当然,有些方法介于这2种方法之间.通常全局模板是指直接从PDB数据库[63]获取的实验测定的蛋白质三维(骨架)结构,而无模板方法是指没有采用全局模板的方法.
TBM方法[14,29~38]大致步骤如图3所示,通常可分为以下几步.第一步,通过数据库检索,得到目标蛋白质的一组同源性序列(MSA),并根据MSA获得1个或多个折叠结构模板.第二步,比对目标序列和模板对应序列,两序列一致的片段直接使用模板的对应折叠结构.第三步,对于目标序列与模板对应序列不一致的区域,采用碎片组装、优化算法或是数据库方法等单独预测.当然通常最后还会用诸如分子动力学的优化方法进行模型的精细化(modelrefinement),以优化全局结构[13,14].历史上,TBM方法[14]可以细分成comparativemodelling(CM)和threading2种方法[29~38].其中在CM中,模板与目标序列的同源性较近.
Fig.3Illustrationofbasicstrategiesoftemplate-basedmethod(TBM)andtemplate-freemethod(FM).
无模板方法[72~77]的流程见图3右半侧,FM从蛋白质数据库中依MSA比对结果找到一些片段的结构并将其放入片段库中,然后找到评分较高的片段结构拼成初始结构[72~74],接着采用FM里非常重要的片段组装(fragmentassembly)方法[72],大致冻结片段的结构并以片段结构为单元来演化全局结构,比如可根据粗粒化的势能函数用梯度下降化进行能量优化.
人们发现在蛋白质变异过程中经常出现关联变异(correlatedmutation)的现象:一条蛋白质链内若发生变异,总是2个氨基酸成对地变异;因为演化压力会迫使蛋白质维持一致构型,原本接触的氨基酸对在变异过程中继续保持接触,可以避免其形状发生剧烈变化.因此,这就使得残基接触对(inter-residuecontactmap)的信息极为重要[56~61].
早期有许多传统方法致力于预测残基接触对.处理该问题的早期算法,倾向于以一次一对的形式、孤立地预测每个接触对是否可能.由于忽视了蛋白质包含的全局信息:一个残基对是否接触受到序列中其他残基的影响,早期算法陷入了困境,预测效果糟糕.而之后研究者提出了充分利用全局信息的预测方法,例如基于Markov随机场模型MRF的directcouplingmethod(DCA)[58~61],在残基接触预测上获得了突破性的成就.
深度神经网络在预测残基接触对问题上,也表现出了异常优异的性能,有时甚至还直接被用于预测键角等信息.这些预测特征均可作为约束,辅助指导无模板方法.
比如,RaptorX-Contact深度学习模型[39]将ContactMap的预测当成图片分割任务来对待,RaptorX-Contact所采用的方法也被其他方法,如ResPRE[54]所采纳.ResPRE[54]采用了图片识别领域非常著名的残差网络(Resnet)模块[70],残差网络的重要思想是不断地将网络前面的信息直接复制到网络后面.
而AlphaFold1[15]又将ContactMap拓展成距离直方图(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.4SketchoftheAlphaFold2model.DetaileddescriptionisreferredtoRef.[12].
而另一输入通道中,主要输入与MSA相对应的序列的结构残基对距离信息以及扭转角的信息.在具体输入时,AF2将距离对长度划分成64个离散块(64bins),并将其转化为概率的形式,故对应数组形状为(s,r,r,64),取值为0~1.注意配对表征中,只包含了MSA除输入序列之外的某个序列自己结构信息,不同序列之间并没有进行信息的关联.
然后再将MSA表征与配对表征输入一个称为Evoformer的模块,此模块主要将MSA的信息(同源性差异)与结构信息整合起来,最后得到输入序列的MSA表征与输入序列的配对表征.此时,输入序列的配对表征同时将演化信息与其他模板结构信息有机地融合在了一起.Evoformer主要利用了自注意力机制来实现上述信息整合.
而下一个结构模块structuremodule主要的功能是将Evoformer预测的配对表征展开成三维空间结构,同时亦承担一定的预测调整功能.此模块的结构大致如图5所示.一条蛋白质骨架结构可想象成一系列三角形的叠加,三角形的中心相当各个残基α碳的坐标,三角形平面本身代表N-α-C-C构成的三角形.这样,此骨架可由2个数组表示,数组形状分别为(r,3×3)和(r,3),分别表示每个三角形取向与位置.
Fig.5Illustrationofhowthepairinginformationistransformedintothe3DstructureusingneuralnetworksinAlphaFold2[12].
初始时,假设所有氨基酸都在原点,然后将此初始骨架与配对表征输入结构模块,由于配对表征存有距离对及取向信息,故可通过一个称为不变点注意力神经网络模块将其初步还原成展开的骨架结构,紧接着再加入侧链原子从而得到全原子的三维结构.
如图4所示,最后再将中间输出的MSA信息、配对信息和3D结构信息重新叠加输入到Evoformer,如此反复迭代3次,最终到预测结果.
因为PDB中只有大约1×105多个的序列有对应的三维结构数据.而在bigfantasticdatabase(BFD)蛋白质序列数据有多达2,204,359,010个序列,虽然这些序列并不一定有对应的三维结构信息(无标签),但self-distillationdataset的训练技巧可以将这些无答案的题目作为作业进行训练,自己提高预测准确度,AlphaFold2用此扩大训练集并进一步提高了预测准确度.
后来有诸多研究团队对AlphaFold2进行了拓展与提升.例如:Baker团队[78]的RoseTTAFold发展了三通路神经网络(three-trackneuralnetwork),对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.6Architectureofthestrongly-correlatedneuralnetwork(SCN)whererindicatesnumberofresiduesandcindicatesnumberoffeaturesorchannels.
该研究发现与传统向量表征方法相比,强关联网络极大提升了预测准确性,提高了约20个百分点.
最近十多年深度学习在蛋白质结构预测中取得了巨大了成就,它的杰出代表AlphaFold2[12]几乎解决半个世纪前提出的蛋白质结构预测难题[12];由于其预测结果达到了中低分辨率的实验精度,几乎等于说AlphaFold2的预测可以直接代替有些蛋白质结构分析实验,而对于通常200多个氨基酸组成的蛋白质,AlphaFold2通常在普通GPU上只需几分钟便能得到其结构,这对于以后的生物制药等领域将有巨大影响.
而另一方面,高分子材料基因组计划仍然在进行中.因为普通高分子的组成不像蛋白质序列那样,有确定的组成单元以及较为单一明确的目标,因此难度更大.但深度学习在蛋白质结构预测中的成功经验仍然对高分子材料基因组计划有一定的启发:
首先,它有一个标准化的结构数据库PDB.高分子材料基因组计划或许也需要构建类似的数据库,难点在于制定统一的数据标准.即如何准确、完整、简洁地表征高分子链,加工条件及性能.
其次,蛋白质结构预测有一个权威的CASP竞赛,CASP极大地推进了结构预测算法的演进.在材料基因组计划中可参照CASP,建立相应的标准化竞赛.
再次,AlphaFold2充分利用了当前深度学习领域的各种先进算法,并不拘泥于某种特定算法.这启发我们解决问题时需要以问题为导向,而非以方法为导向.
最后,AlphaFold2中将Distogram信息用神经网络转化成分子结构坐标的方法可推广至其他结构预测的问题中,当然也可用于高分子的结构预测.