1、引言

蛋白质在生命体的基础生物学活动中扮演着重要角色，是生物体内最重要的分子之一，它们在细胞内执行许多关键的功能，包括催化化学反应、传递信号、提供结构支持和运输分子等[1,2]。通过计算机技术预测蛋白质三维结构可为药物设计、疾病诊断和精准治疗提供重要的理论基础[3]。

模板对于蛋白质结构预测是非常重要的。传统的蛋白质结构预测包括模板建模(Template-based Modeling)和从头预测(Ab Initio Modeling)两种方法[4,5]。模板建模方法是利用已知的蛋白质结构作为模板来预测目标蛋白质的结构[6]。通过将目标蛋白质序列与已知结构的库进行比对，找到相似的模板结构，并将模板的结构信息映射到目标蛋白质上，生成初始的结构模型。常见的模板建模方法包括SWISSMODEL[7],RosettaCM[8]和MODELLER[9]。这些方法都依赖于已知的蛋白质结构作为模板，通过比对、拼接、约束等方式来预测目标蛋白质的结构。蛋白质结构从头预测法是在没有合适模板的情况下，从目标蛋白质的序列出发，通过搜索蛋白质构象空间中的最低自由能状态，来寻找最稳定的结构[10]。常见的蛋白质结构从头预测方法包括D-I-TASSER[6],Rosetta[11]和MMpred[12]。这些方法虽然没有显式地使用模板，但它们都使用了基于模板结构构建的片段信息[13]。这些片段信息可以提供有关蛋白质的结构特征和构象约束。

除了传统的建模方法，基于机器学习的蛋白质预测技术也十分依赖于模板结构[14]。基于机器学习的蛋白质预测方法通常通过训练模型来学习蛋白质序列和结构之间的对应关系[6,15]。这些方法可以利用已知的蛋白质结构作为训练数据，提取结构特征和序列特征来建立预测模型[15]。此外，还可以结合传统的模板建模方法，通过将模板建模和机器学习技术相结合，综合利用多种信息源来进行结构预测，包括AlphaFold2[16],RoseTTAFold[17]和trRosettaX[18]。这些方法将模板结构作为输入特征之一，与其他结构特征和序列特征一起输入机器学习模型中进行训练和预测，提高预测的准确性和泛化能力[19]。因此，开发一种高效的方法来识别高质量的远程同源模板对蛋白质结构预测是至关重要的。

总体来讲，现有的模板检测方法包括3类，分别是序列-序列比对、谱-谱比对和穿线法。序列-序列比对的模板检测方法包括BLAST系列[20]。BLAST是将目标蛋白质序列与数据库中的已知结构蛋白质序列进行比对搜索，根据相似性得分和期望值来评估比对结果。PSI-BLAST[21]是BLAST的一个改进版本，它利用迭代比对的方式增强序列相似性搜索的灵敏度。谱-谱比对的模板检测方法包括HHsearch[22],MUSTER[23]和 SPARKS-X[24]等。这些方法使用了序列谱信息、二级结构、溶剂可及性和扭转角来构建评分函数，相比序列-序列比对，提高了模板搜索的精度和覆盖率。穿线法包括EigenTHREADER[25]和CEthreader[26]等。这类方法首先预测目标序列的残基接触图并转化为特征向量表示，然后通过特征向量分解和动态规划算法的结合，可以在结构空间中寻找最佳的结构对齐。除此之外，元穿线法LOMETS系列[27,28,29]集成了各类的模板检测方法，利用多种模板检测方法的优势，提高了蛋白质结构预测的准确性和覆盖范围。虽然这些方法在一定程度上提升了模板检测的精度，但仍然有很大的提升空间。

本文设计了一种基于特征向量自适应提取比对的远程同源模板检测算法(Adaptively Selected Eigenvector alignment, ASEalign)。其采用多特征信息融合的深度学习技术预测蛋白质接触图，并设计了多维度特征打分函数，通过自适应地提取接触图矩阵中的特征值和特征向量进行模板比对，提升了模板检测的效率和准确性。通过将检测出的模板输入AlphaFold2中进行结构建模，提升了AlphaFold2的模型精度，进一步验证了模板结构对蛋白质建模的重要性。

2、方法

远程同源模板检测算法ASEalign的流程图如图1所示。从序列出发，通过HHblits[30]搜索UniRef30[31]和BFD库生成MSAs并从中提取一维和二维特征，通过水平条带化将一维特征和二维特征进行组合，得到一个L×L×490的张量，然后将特征张量输入自注意力机制模块和卷积残差模块中预测出残基接触图。基于接触图、二级结构、序列谱谱比对和溶剂可及性设计的打分函数，自适应选取接触图矩阵中的特征值和特征向量与PAcluter80[32]模板库进行比对，检测出最终的模板结构。最后使用AlphaFold2预测器进行模板结构建模。

图1 ASEalign流程图

2.1 蛋白质接触图预测

2.1.1 训练集和测试集的构建

PDB库1)收集了所有已被生物实验方法测定的蛋白质结构，提供了丰富的模板结构信息[33,34]。本文使用PDB数据库构建了用于蛋白质接触图预测的训练集。首先从PDB数据库中挑选由X射线衍射解析出的蛋白质，然后筛选出序列长度在40～500之间且分辨率小于2.5Å的蛋白质作为候选数据。使用CD-HIT工具[35]对收集到的蛋白质进行聚类，以30%的序列相似度为阈值去除冗余，最终得到13 980条蛋白质序列作为训练集，其中95%用于训练，5%用于验证。

为了更加客观地评价方法的性能，本文从SCOPe 2.07[36]数据库中构建了测试集[37]。首先通过CD-HIT工具以30%的序列相似度为阈值去除冗余，提取到11 198条蛋白质序列；然后选取长度在50～500之间且分辨率小于1.5Å的蛋白质作为候选数据；再一次使用CD-HIT工具以30%的序列相似度对训练集进行去冗余处理，最终得到135个序列作为测试集。

2.1.2 特征提取

首先为训练集中的每一条蛋白质序列搜索MSA。利用HHblits, 将其E-value的阈值设置为0.001,序列覆盖率设置为至少50%,在UniRef30和BFD为训练集的蛋白质序列生成MSAs, 并且从该目标序列以及MSAs中提取蛋白序列的一维特征和二维特征。一维特征包括序列频率谱、溶剂可及性信息、二级结构信息；二维特征包括协方差矩阵、CCMpred耦合分数[38]、残基对接触势能、去除背景噪声的互信息。

序列频率谱：表示残基在MSAs中出现的频率。不同的频率表征了氨基酸在蛋白质序列中的位置和特性。

溶剂可及性信息：指氨基酸在溶液中受到水分子包围的程度。溶剂可及性反映了氨基酸在溶剂中亲疏水性所引起的相互作用力。

二级结构信息：蛋白质的折叠过程中形成的局部稳定的α螺旋、β折叠等结构单元。

协方差矩阵：是通过MSA中的残基频率计算得到的，描述了蛋白质序列中任意两列残基之间的边缘分布和联合分布之间的相关性。其对角线元素表示残基的方差，非对角线元素表示不同残基之间的协方差，反映了结构中的相互作用。

CCMpred耦合分数：CCMpred描述残基间的非线性关系，通过马尔可夫随机场(Markov Random Field, MRF)学习MSA数据的生成模型，根据残基对接触势能的Frobenius范数应用平均乘积修正来消除传递相互作用。

残基对接触势能：不同氨基酸残基之间的接触势能，主要包括范德华力、氢键作用、静电吸引力等相互作用力引起的势能变化。

去除噪声的互信息：互信息用来表征MSA中固定位置残基在共进化过程中的共变程度。为了消除进化压力对氨基酸出现频率的影响，利用氨基酸背景频率修正互信息中的边缘频率分布，引入去除噪声的互信息。

2.1.3 网络搭建

采用4层注意力机制以及128个卷积残差块的网络模型。注意力机制可以对不同的残基特征进行加权学习和融合，使得模型能够更加关注那些对预测接触关系最为重要的特征，捕捉到不同残基的长程依赖关系。卷积神经网络通过卷积以滑动一个卷积核的形式来提取蛋白质特征，利用池化操作对特征图进行降维，再利用全连接层将卷积层和池化层提取的特征进行组合，最后引入非线性因素拟合更复杂的数据。多头注意力机制和卷积神经网络能够自动学习和提取蛋白质序列中的残基间进化关系，并输出其映射到的接触图。两者相互协作，更加全面地描述蛋白质的特征，提高了模型的泛化能力。

首先把提取的一维特征和二维特征进行组合，得到一个L×L×490的特征张量。通过二维实例归一化层对特征数据进行规范化处理后，得到L×L×128的特征张量，输入两层自注意力机制模块中，使用4头注意力同时关注多个空间特征，防止单一注意力机制将信息集中于自身。然后将得到的L×L×128张量输入卷积模块。卷积模块包含128个残差块，每个残差块包含一个二维卷积层(卷积核为3*3)、一个批量归一化(BP)层、一个指数线性单元(ELU)激活层、一个dropout层(dropout rate 20%)和一个二维卷积层(卷积核在1,2,4,8,16核之间交替膨胀)。

2.1.4 模型训练

网络模型是用于预测目标蛋白质残基之间的接触图，残基之间的是否接触定义为残基Cβ原子(甘氨酸为Cα原子)之间的距离是否小于等于8Å。在训练阶段，从训练集中选择95%的蛋白质用于训练，5%的蛋白质用于验证。为了防止训练过程中模型会形成记忆，在每一轮训练之前都会打乱蛋白质的顺序。训练时在序列距离大于4的残差对上计算损失，并使用预测和真实接触点之间的二元交叉熵作为损失函数，该网络模型的训练共迭代了50个epoch。

2.2 自适应接触矩阵特征分解

自适应接触矩阵特征分解是一种自适应选取接触矩阵的特征值和特征向量进行模板比对的算法。当矩阵的特征值较大时，特征空间所包含的信息量也较大。通过特征向量相关性最大的几个特征向量来近似接触矩阵，减少了比对次数。自适应特征向量比对方法能够根据特征值、特征向量灵活调整其贡献率，在保证算法精度的同时，提高了算法的效率。与固定特征值数量的模板比对算法相比，自适应选取特征向量的模板比对方法能保留较长序列蛋白质更多的特征值，具有较大的优势。因此，用两组少量特征向量的全局对齐计算两个接触图的重叠程度，能有效地提高匹配的准确度和鲁棒性。

接触矩阵M是对角线均为0的实对称矩阵，有L个特征值和相对应的特征向量，根据谱图理论可知：

M=VΔV-1=VΔVT (1)

其中，λi表示第i个特征值，i∈{1,2,…,L}以及

表示中对应的特征向量。将特征值按递减顺序排序为λ1≥λ2≥λ3≥…≥λi≥…≥λL。根据式(2)、式(3),对接触矩阵全部特征值总和的贡献比例进行如下计算：

其中，C表示λi的贡献率。将贡献率累加计算，为使得接触矩阵贡献率之和刚好大于阈值的第一个特征值的索引值，k∈{1,2,…,L}。其中，C是0.5的阈值。利用筛选出的k个特征值及其对应特征向量重塑接触矩阵得到接触矩阵的近似表达如式(4)所示：

则接触矩阵M残基i和j之间的接触信息可以近似为式(5):

其中，[vi,1 … vi,k… 0]表示自适应选出i位置的特征向量，i∈{1,2,…,L},[vj,1 … vj,k … 0]表示自适应选出位置j的特征向量，j∈{1,2,…,L}。

2.3 模板检测打分函数设计

不同于传统的仅利用序列信息进行模板比对，本文根据多维度特征设计了更加精准的打分函数，通过多特征融合的评分函数来提高模板库筛选的准确性。使用动态规划算法进行序列-模板对齐打分。

首先，将目标蛋白p1和模板蛋白p2中的残基分别表示为R1={1,…,n};R2={1,…,m}。接触图分别是MP1∈{0,1}n×n和MP2∈{0,1}m×m。目标蛋白的第i个残基与模板蛋白的第j个残基接触匹配得分的计算式如(6)所示：

其中，R=3,E1(i,j),E2(i,j)和E3(i,j)分别表示序列谱谱比对得分、二级结构得分和溶剂可及性得分。详细介绍及参数可参见文献[39]。 Ec(i,j)表示目标蛋白的第i个残基与模板蛋白的第j个残基的通过接触匹配的得分，计算式如式(7)和式(8)所示：

Ec(i,j)=Econ(i,j)+Egap(i,j) (7)

分别表示p1和p2中的第i和第j个残基与其他残基形成接触的特征向量；k和t分别表示自适应选取后特征向量的个数，1≤k≤n, 1≤t≤m; Egap(i, j)是比对过程中加入Gap的惩罚分数。

3、结果分析

3.1 模板检测精度比较与分析

为了检验ASEalign方法的预测性能，本文在135个蛋白质上进行了模板检测性能的测试，并与基于隐马尔可夫的谱-谱比对方法HHsearch以及基于接触图比对的EigenThrea-der进行了结果对比。使用TM-score评估了模板的精度，如图2所示。它考虑了两个结构的全局拓扑特征。TM-score的取值范围在0～1之间，值越接近1表示两个结构越相似，值越接近0表示两个结构越不相似。ASEalign在135个测试蛋白上检测模板的TM-score均值为0.695,比HHsearch(0.623)高出11.5%,比EigenThreader(0.636)高出9.2%。当TM-score ≥ 0.5时，模板结构与天然蛋白的拓扑结构非常相似。统计两个算法检测的模板的TM-score ≥ 0.5的数量，ASEalign有110个，占总测试集的85%。这表明ASEalign的模板检测性能有了显著的提升。在135个测试蛋白质中，ASEalign有85个蛋白质的模板检测结果好于HHsearch。其余50个模板检测结果比HHsearch差，主要是因为ASEalign只搜索了PAcluster80的质心结构，这提升了模板检测的速度，但同时也损失了一部分模板的精度。

图2 ASEalign 与 HHsearch 和EigenThreader检测模板的TM-score 箱型图

图3 ASEalign 和 HHsearch 在肽基-tRNA水解酶上检测出的模板与天然结构的比较

本文对ASEalign和HHsearch检测的模板进行了详细的比较与分析。肽基-tRNA水解酶是一种细菌酶，可裂解肽基-tRNA或N-酰基-氨酰-tRNA以产生游离肽或N-酰基-氨基酸和tRNA。肽基-tRNA水解酶是一个单一的α/β球状结构，具有7个β链，形成一个扭曲的中心β折叠，被6个螺旋包围。对于HHsearch, 它检测出的模板精度为0.398,仅仅包含了中心的β折叠和一个α螺旋，没有识别出具有完整结构区域的模板。而ASEalign检测出的模板精度为0.721,该模板基本覆盖了整个目标蛋白，结果显著好于HHsearch。这是因为该目标蛋白的MSAs数量仅有653条，HHsearch基于少量的MSAs提取的隐马尔可夫谱信息是有限的，这降低了HHsearch的模板精度。而ASEalign使用同样数量的MSA信息，通过机器学习预测出了准确的接触图，通过提取特征向量进行比对检测出了精度更高的模板。这表明ASEalign基于接触图比对的模板检测能够比HHsearch基于谱-谱比对的模板检测获得更精确的模板结构。

为了进一步检验ASEalign远程同源模板的检测性能，本文在一周的CAMEO数据集(2023/08/19)上进行了远程模板(去除了大于等于30%序列相似的同源模板)检测并与HHpred进行了比较，如图4所示，共包含了16个蛋白质。ASEalign检测的模板的平均TM-score为0.796,比HHpred的平均TM-score高出5%。该结果再次表明ASEalign对远程同源模板的检测性能优异。

图4 ASEalign与HHpred在CAMEO数据集上检测的远程模板的比较

3.2 自适应选取特征向量的结果分析

为了提高模板检测的准确性并加快模板检测速度，采用了自适应选取特征向量的策略进行模板比对，而不是固定数目的特征向量。特征向量的数目是根据特征值贡献率来选择的。特征值贡献率是接触矩阵每个特征值与所有特征值总和的比值。本文分别对0.3～1.0的特征值贡献率进行了实验，如图5所示。

图5 特征值贡献率与特征向量个数、运行时间和TM-score的关系

随着贡献率的逐步提高，特征向量的个数和模板精度(TM-score)都逐步增加。这表明越多的特征向量所包含的信息越丰富，将有助于提升模板检测的精度。然而，当贡献率增加时，模板检测的时间也随之增加。本文选择了特征值贡献率为0.8作为最终比对方案，既在0.8特征值贡献率阈值下，采用自适应选取前80%的特征值及其对应的特征向量用于模板比对。在这个阈值下，可以保证特征向量比对数量在趋势激增之前，既保留了包含主要接触信息的特征值，避免浪费非必要的特征比对，同时减少了模板检测所运行的时间。对于长度不同的蛋白自适应地选取不同特征数量的信息进行比对，提高了模板检测的精度和效率。

为了检验自适应选择特征向量的有效性，本文与CEthreader做了进一步的比较，如表1所列。CEthreader采用固定数量为7的特征向量进行模板比对；ASEalign采用自适应选取特征向量比对的策略，在所有目标蛋白上选取的特征向量平均个数为49。结果显示，ASEalign的平均TM-score为0.695,比CEthreader高出了8.4%。这是因为在一些较长的蛋白中，固定的特征值数量可能会导致搜索模板时丢失一些额外的信息，从而降低模板比对的精度。相比之下，ASEalign自适应选取的特征值基本覆盖了蛋白质残基间的接触信息，因此在保证搜索效率的前提下提升了模板检测的精度。

表1 ASEalign和CEthreader的模板检测结果

3.3 模板增强的AlphaFold2建模

为了进一步验证ASEalign检测的模板的有效性，本文通过模板建模与AlphaFold2进行了比较。AlphaFold2使用HHsearch搜索的前4个模板进行建模。本文用ASEalign检测的模板替换掉HHsearch的模板，并使用AlphaFold2进行了建模，并将方法命名为ASEalign-AF2,结果如表2所列。AlphaFold2预测的模型平均TM-score为0.834,通过ASEalign-AF2模板增强的模板建模精度提升到了0.839,较初始模型提升了0.6%。在所有目标蛋白质中，其中有76个的结果比初始模型效果更好，占比58%。在RMSD的比较上，经过模板增强后的AlphaFold2模型的平均RMSD值从5.42Å下降到了5.14Å。这表明通过模板增强的建模确实可以提高AlphaFold2的预测准确性。

表2 ASEalign-AF2和AlphaFold2的建模结果

以两域蛋白质6J93_A为例，本文对ASEalign-AF2和AlphaFold2的建模结果进行了分析，如图6所示。ASEalign-AF2的建模精度为0.893,AlphaFold2的建模精度为0.836。ASEalign-AF2的模型精度显著高于AlphaFold2,这主要得益于ASEalign提供了准确的模板信息。ASEalign检测的模板结构的TM-score为0.788,而HHsearch检测的模板结构的TM-score为0.431。从图6可以看出，HHsearch检测的模板在单域上的精度较好，而域间的方向是不正确的。ASEalign检测的模板提供了正确的域方向，因此得到了较高精度的模型。这表明ASEalign可以提供准确的模板结构信息用于模板建模。

图6 在示例蛋白质6J93＿A上ASEalign和Alpha Fld2o的模板和建模结果

结束语蛋白质的结构对于其功能和相互作用至关重要，因此准确地预测蛋白质的结构是理解其功能和设计新药物的关键一步。模板结构为蛋白质的预测提供了一个框架，可以指导目标蛋白质的结构建模。尤其是对于没有同源序列的孤儿蛋白，它的结构预测将十分依赖于远程同源模板。因此，检测出合适的模板结构对于预测高精度的蛋白质结构非常关键。

本文设计了一种基于蛋白质接触图自适应选取特征向量比对的远程同源模板检测算法ASEalign。从序列出发，通过HHblits搜索MSAs并从中提取一维和二维特征，通过水平条带化将一维特征和二维特征进行组合，得到一个L×L×490的张量，然后将特征张量输入自注意力机制模块和卷积残差模块中预测出残基接触图。基于设计的多维度特征打分函数，自适应地选取接触图矩阵中的特征值和特征向量与PAcluter80模板库进行比对，检测出最终的模板结构并使用AlphaFold2预测器进行模板增强的结构建模。在135个蛋白的测试集上的结果表明，ASEalign相比主流的模板检测算法HHsearch精度提升了11.5%,并通过模板增强的蛋白质结构建模提升了 AlphaFold2的单体模型精度。这表明ASEalign检测远程同源模板的性能优于目前主流的模板检测算法。

尽管AI技术的发展给单域蛋白质结构的预测带来了巨大的进步，但多域蛋白质建模以及复合物的组装仍然存在很大的挑战，这些都离不开模板提供的结构信息[40]。除此之外，模板检测在蛋白质功能注释、药物设计以及生物工程和蛋白质工程等领域都具有重要的用途[19,41]。它们提供了有关蛋白质结构和功能的宝贵信息，为科学研究和应用开发提供了基础和指导。

基金资助:国家自然科学基金(62173304); 国家重点研发计划(2019YFE0126100)~~;

文章来源:梁方,徐旭瑶,赵凯龙,等.远程模板检测算法及其在蛋白质结构预测中的应用[J].计算机科学,2024,51(S1):179-185.