考虑转移熵和空间近邻进化信息的基于集成模型的蛋白质变构位点预测方法与流程

1.本发明属于蛋白质功能位点预测技术领域，是一种基于蛋白质三维结构信息的变构位点预测方法。


背景技术：

2.变构是蛋白质发挥功能的重要调控手段，通常是指调节分子(变构剂)结合于蛋白质活性部位(正构位点)以外的其它部位(变构位点)，引起活性部位结构或动力学性质的改变，从而影响蛋白质功能的现象
[1-3]
。比较正构与变构位点的空间和拓扑结构显示变构位点的结构更具多样性，这使得变构调节剂具有更高选择性和更少毒副作用，这些优势使得变构剂设计在药物研发中具有更重要的意义
[4-6]
。变构位点的发现是变构药物设计的前提，如何准确高效地识别蛋白质变构位点是结构生物学家和药物设计者普遍关心的问题。
[0003]
目前，实验上获得变构分子/变构位点多是由高通量筛选实验偶然所得
[7-8]
，这极大阻碍了变构药物的设计与开发
[9]
。发展理论方法来准确预测蛋白质变构位点是亟待解决的问题。当前，许多理论方法被提出用以研究蛋白质变构及识别变构位点。其中，基于分子动力学的经典方法有双态模型
[10]
、allomod
[11]
和spacer
[12]
等方法。这些方法基于物理学思想，需要对构象系综进行充分采样来发现可能的变构位点，通常需要耗费较长的时间。对此，人们又提出了基于正则模分析和粗粒化弹性网络模型的方法来识别蛋白质变构位点，包括pars
[13]
和corrsite
[14]
等方法。除了动力学的方法，研究人员还提出了基于进化分析的方法
[15-16]
和基于结构网络的方法
[17-18]
，但预测准确性都不够理想。对此，人们开始试图建立可以考虑多因素的机器学习方法来进行预测。2013年，huang等人首次将机器学习方法用于变构位点预测，提出了allosite方法
[19]
。该方法使用fpocket
[20]
查找蛋白质表面潜在的变构口袋，并将其理化特性作为口袋的描述符来训练模型。2017年，song等人利用口袋描述符和蛋白质与虚拟配体结合时柔性变化的显著性作为特征，开发了allositepro方法
[21]
。2022年，xiao等人开发了passer 2.0网络服务器用于变构口袋的识别
[22]
，该方法利用autokeras
[23]
和autogluon
[24]
算法来解决特征工程、模型选择和超参调优等耗时费力的问题。
[0004]
最新研究表明，蛋白质变构活性在本质上是熵的效应，它依赖于由残基间协同涨落引起的信息从变构位点向正构位点的传递
[25-26]
。近来，由schreiber
[27]
提出的转移熵概念，由于其与信息流和变构交流有关，得到了广泛关注。转移熵中引入了残基间时间延迟运动相关性信息，因此它可以呈现并量化残基间驱动与被驱动的关系，从而揭示哪些位点可以被操控来影响蛋白质正构位点的活性，这些位点很可能成为蛋白质的变构位点。另外，蛋白质功能位点残基，如配体结合和催化位点，往往不是独立发挥作用，而是彼此协同发挥功能
[28-30]
，蛋白质变构口袋残基也一样。这种协同性必然会体现在残基进化上，从而呈现出一定的空间协同进化特点。最近，我们提出了一种新的进化信息编码方法，即基于空间近邻的位置特异性打分矩阵(spacial neighbor-based position-specific scoring matrix,
residues at protein-rnainterfaces and its application in binding site identification[j].bmc bioinformatics.2020,21(1):57.
[0034][0035]
[29]guharoy m,chakrabarti p.conserved residue clusters at protein-protein interfaces and their usein binding site identification[j].bmc bioinformatics.2010,11:286.
[0036]
[30]ahmad s,keskin o,sarai a,et al.protein-dna interactions:structural,thermodynamic andclustering patterns of conserved residues in dna-binding proteins[j].nucleic acids res.2008,36(18):5922-5932.
[0037]
[31]liu y,gong w,yang z,et al.snb-pssm:a spatial neighbor-based pssm used forprotein-rna binding site prediction[j].j mol recognit.2021,34(6):e2887.
[0038]
[32]liu y,gong w,zhao y,et al.aprbind:protein-rna interface prediction by combiningsequence and i-tasser model-based structural features learned with convolutional neuralnetworks[j].bioinformatics.2021,37(7):937-942.


技术实现要素：

[0039]
本发明的目的是在已知变构蛋白质结构的情况下，预测其变构位点，为蛋白质变构位点识别和变构药物设计提供重要信息。方法的特点在于：
①
查找蛋白质表面潜在的变构口袋，并提取口袋特征。除口袋的理化性质外，本发明还提取了转移熵和残基协同进化信息。结果显示这些新加入的特征对变构位点预测有好的贡献，在特征选择中被选入最优特征组合。
②
本发明提出了基于多份训练子集的集成模型构建策略，即拆分大类样本、产生多份训练子集、选择适合多份训练子集的最优特征组合、构建基于多个子模型的集成模型来进行变构位点预测。该策略由于利用了训练集中所有的大类样本，因此它除了有利于模型的泛化能力提高外，还附带解决了训练集中变构位点正负样本极度不平衡的问题。
③
本发明提出了基于多份训练子集的集成模型构建策略，并首次将转移熵和空间近邻进化信息用于蛋白质变构位点预测，其中空间近邻进化信息是我们之前开发的(jmol recognit.2021,34(6):e2887)，可以很好地考虑氨基酸残基空间近邻的协同进化性。
[0040]
本发明考虑转移熵和空间近邻进化信息的基于集成模型的蛋白质变构位点预测方法，是一种基于转移熵和空间近邻进化信息的集成分类模型来对蛋白质变构位点进行预测的方法(ensemble-based machine learning approach with transfer entropy and snb-pssm considered for protein allosteric site prediction,allosests)，包括四个步骤：一是查找蛋白质表面潜在的变构口袋，二是提取口袋特征，三是产生多份训练子集并筛选获得最优特征组合，四是构建基于多个子模型的集成模型以预测蛋白质变构位点(如附图1)。
[0041]
建立模型前，首先构建了训练集和测试集：
[0042]
收集来自蛋白质数据库(protein data bank,pdb，https://www.rcsb.org/)、变构数据库(asd,https://mdl.shsmu.edu.cn/asd/)、变构基准数据集(asbench,https://
mdl.shsmu.edu.cn/asbench/)和song等人在2017年提出allositepro方法的工作(j chem inf model.2017,57(9):2358-2363)中构建的测试集中的变构蛋白质。对这些变构蛋白质，按照序列一致性小于30％、结构分辨率优于的标准筛选出179个变构蛋白质作为数据集，其中来自asbench数据集的127个蛋白质作为训练集，来自allositepro工作的24个蛋白质作为测试集1(test set 1)，来自从多个文献中收集整理的28个蛋白质作为测试集2(test set 2)；表1-3分别列出了上述三个数据集中的变构蛋白质的信息，包含pdb id、链号、残基数和口袋样本数；
[0043]
表1训练集中的127个变构蛋白质及其口袋样本数
[0044]
[0045]
[0046]
[0047]
[0048][0049]
a/b:变构口袋数/口袋总数。
[0050]
表2测试集1中的24个变构蛋白质及其口袋样本数
[0051]
[0052][0053]
a/b:变构口袋数/口袋总数
[0054]
表3测试集2中的28个变构蛋白质及其口袋样本数
[0055]
[0056][0057]
a/b:变构口袋数/口袋总数
[0058]
下面对蛋白质变构位点预测方法的四个步骤进行描述。
[0059]
步骤1：查找蛋白质表面潜在的变构口袋
[0060]
针对选定的变构蛋白质，本发明保留带有变构残基的功能链，去除水分子和配体，将其用于查找蛋白质表面潜在的变构口袋，这些变构口袋被作为样本数据；样本数据中与变构调节剂原子接触最多的口袋被标记为阳性样本，其余被标记为阴性样本；对于含有127
个蛋白质的训练数据集，查找到2245个口袋，其中有147个阳性样本和2098个阴性样本；对于两个测试集，即包含24个蛋白质的测试集1(test set 1)和包含28个蛋白质的测试集2(test set 2)，分别查找到24和28个阳性样本，295和347个阴性样本；
[0061]
这里，使用高效开源的fpocket(bmc bioimformatics 2009,10:168)工具查找蛋白质表面潜在的变构口袋，采用默认的参数设置；对于查找到的蛋白质表面口袋，fpocket提供了19个物理化学特征，如表4所列。这些口袋如前所述被分成训练集和测试集用于预测模型的训练和测试；
[0062]
表4fpocket计算的口袋理化性质特征
[0063][0064][0065]
步骤2：提取蛋白质口袋特征
[0066]
对查找到的蛋白质表面潜在变构口袋，除19个理化性质特征(表4)外，还提取了16个表征口袋二级结构信息的特征、20个表征空间协同进化信息的特征，和1个表征动力学转移熵信息的特征。每一个蛋白质表面口袋对应的后三类特征的提取具体描述如下：
[0067]
口袋二级结构信息的特征：是根据其中残基的8种二级结构类型来提取的，包括口
袋中氨基酸残基在每种二级结构类型中出现的数量，和该数量占口袋内全部残基数量的百分比。8种二级结构类型包括：α-helix(h)、β-bridge(b)、β-ladder(e)、3
10-helix(g)、π-helix(i)、turn(t)、bend(s)和不确定类型(m)，二级结构类型使用dssp工具(biopolymers.1983,22(12):2577-2637)计算；
[0068]
空间协同进化信息的特征：我们之前(jmolrecognit.2021,34(6):e2887)提出的基于空间近邻的位置特异性打分矩阵(snb-pssm)被用于提取蛋白质口袋的进化信息。首先，对目标蛋白质序列，用psi-blast工具(nucleicacidsres.1997,25(17):3389-3402)在非冗余蛋白质序列数据库中循环查找其同源序列，以e-value＝0.001为阈值迭代运行三次，生成标准位置特异性打分矩阵(position-specificscoringmatrix,pssm)。pssm给出了同源序列比对的结果，其中包含进化信息；对于含有n个残基的蛋白序列，pssm矩阵的大小为20
×
n，每个残基位置的进化信息被编码在20维的向量中；然后，为了考虑目标残基空间近邻的进化信息，对其原始进化信息进行基于空间近邻的平滑处理，使其进化信息为：与该残基ca原子距离以内的所有残基进化信息的平均值(见附图2的一个例子)。最后，口袋内所有残基进化信息的平均值被作为描述口袋的空间协同进化信息特征，该特征为一个20维的向量。
[0069]
动力学转移熵信息的特征：动态高斯网络模型(dynamicgaussiannetworkmodel,dgnm)(proteins.2017,85:1056-1064)被用于计算转移熵；蛋白质结构是动态的，平衡态下蛋白质在围绕其平均构象涨落运动，t时刻蛋白质构象偏离其平均构象可用一个3n维的矢量δr(t)来表示(用c
α
原子的位置记录一个残基的位置，n为残基数)。根据schreiber的转移熵理论，从残基i的轨迹到残基j的轨迹(时间间隔为τ)的转移熵ti→j(τ)为：在给定残基j轨迹t时刻对应的过去值δrj(t)，知道过去值δri(t)时，未来t τ时刻δrj(t τ)值不确定性的减少量：
[0070]
ti→j(τ)＝s(
△rj
(t τ)|
△rj
(t))-s(
△rj
(t τ)|
△ri
(t),
△rj
(t))(1)
[0071]
这里，s(δrj(t τ)|δrj(t))是给定δrj(t)时，δrj(t τ)的条件熵；s(δrj(t τ)|δri(t),δrj(t))是给定δri(t)和δrj(t)时，δrj(t τ)的条件熵；根据香农熵，等式(1)中的两项为：
[0072]
s(
△rj
(t τ)|
△rj
(t))＝-《lnp(
△rj
(0),
△rj
(τ))》 《lnp(
△rj
(0))》(2)
[0073][0074]
将等式(2)和(3)代入(1)得到：
[0075][0076]
其中，《》表示系综平均，p()表示概率。
[0077]
在gnm模型下(以c
α
原子为节点，截断半径为时间延迟交叉相关和转移熵分别为：
[0078]
[0079][0080]
这里，和分别是第k个运动模式下残基i与j的交叉相关性、残基j的自相关和残基i的自相关，τ0为特征时间(characteristictime)，这里采用τ/τ0＝5的时间延迟。运动模式由gnm下蛋白质体系的基尔霍夫矩阵获得，基尔霍夫矩阵为：
[0081][0082]
其中，为截断半径，r
ij
表示两个节点间的距离。基尔霍夫矩阵的伪逆矩阵为：
[0083][0084]
其中，k为运动模式数，λk为对应第k个运动模式的特征值，μk为对应第k个运动模式的特征向量。
[0085]
从残基i到j的净熵转移(netentropytransfer,net)是从i到j的熵转移与从j到i的熵转移之差：
[0086]
neti→j＝ti→j(τ)-tj→i(τ)(9)
[0087]
从残基i到蛋白质其他残基的净熵转移neti→e是上式对残基j求和。为方便比较，对neti→e进行归一化处理。口袋的净熵转移被定义为口袋中所有残基净熵转移的平均值，该值若大于0，则口袋倾向于发出信号，否则倾向于接收信号；口袋的净熵转移即为口袋动力学转移熵信息的特征。
[0088]
步骤3：产生多份训练子集并筛选获得最优特征组合
[0089]
在变构位点预测中，步骤1中对应的正负样本数量是极不平衡的，可达到1:10。如此高的不平衡将严重影响分类器模型的预测性能，造成对大类样本的过拟合。通常，过采样
或欠采样技术被用于处理不平衡的数据。过采样是通过随机复制或生成小类样本，而欠采样则是通过删除部分大类样本，来达到两类数据的平衡。两种技术分别可能会产生不合理的样本或造成一定程度信息的损失。为此，我们提出了一种基于多份训练子集的集成模型构建策略来克服样本的不平衡问题。首先，将训练集中的大类样本随机分为10个子集(子集中样本数与小类样本数相当)，并将它们分别与小类样本组合成10份训练数据集；接下来，对每一份训练数据集，使用前向特征选择算法(sequential forward selection，sfs)从预先提取的特征(步骤2中的56种特征)中筛选出最优特征子集；这里，我们通过对每一份训练数据集进行5次10-fold交叉验证来获得最优特征子集，然后将10个最优特征子集取并集获得一个特征子集；对该特征子集，使用最大相关最小冗余算法(maximum correlation minimum redundancy,mrmr)进一步筛选获得最终的最优特征组合；mrmr能够保证所选特征与输出变量间的最大相关，同时保证所选特征间的最小冗余；表5列出了最终筛选出的最优特征组合；
[0090]
表5特征选择筛选的最优特征组合
[0091]
[0092][0093]
步骤4：构建基于多个子模型的集成模型以预测蛋白质变构位点
[0094]
针对每一份训练数据集，建立相应的支持向量机(support vector machine,svm)分类子模型，通过5-fold交叉验证来训练模型；然后对10个子模型的预测结果(即样本是变构口袋的概率)进行平均整合，构建集成分类模型，来预测口袋是否为蛋白质的变构口袋(如附图3的一个例子)。输入数据为表5列出的最优特征组合数据，输出为预测的蛋白质变构口袋的位置和概率排名；
[0095]
allosests模型在10份训练集上均使用5-fold交叉验证训练和调参(通过最大化马修斯相关系数mcc)，并在独立测试集1和2上进行测试；模型的评价指标包括敏感性(sensitivity,sen)、特异性(specificity,spe)，精度(precision,pre)和马修斯相关系数(matthews correlation coefficient,mcc)，其公式如下：
[0096][0097][0098][0099][0100]
其中，tp(true positive)、fp(false positive)，tn(true negative)和fn(false negative)分别为真阳性、假阳性，真阴性和假阴性样本数量，它们通过比较每个口袋的预测情况与真实情况来确定。对于模型的测试，我们将测试蛋白质直接输入到集成模型。
[0101]
本发明提出了一种考虑转移熵和空间近邻进化信息的基于集成模型的蛋白质变构位点预测方法allosests。该方法能够比较快速地预测蛋白质变构口袋的位置并进行排名，在预测性能上比现有的其他预测方法具有更好的评价指标。与其他预测方法相比，本发明不仅考虑了蛋白质口袋的理化性质和结构特征，还考虑了动力学转移熵信息和残基协同进化信息，是一个基于多个子模型的集成分类器。allosests在蛋白质变构口袋预测中有较好性能，可用于蛋白质变构药物靶点识别和变构药物设计的研究。
附图说明
[0102]
图1为考虑转移熵和空间近邻进化信息的基于集成模型的蛋白质变构位点预测方法流程图；
[0103]
图2为snb-pssm方法中目标残基进化信息编码的过程；(a)目标残基5的空间近邻残基(3、4、6、13和14)，它们的cα原子在残基5的cα原子的范围内。(b)残基5的进化信息被定义为其空间近邻残基3、4、6、13和14在标准pssm中进化信息(20维的向量)的平均值。
[0104]
图3为本发明所采用的基于多个子模型的集成分类器预测框架图；
[0105]
图4为allosests预测的变构口袋位置及口袋排名。
具体实施方式
[0106]
下面结合实施例对本发明做进一步的说明，但本发明并不限于以下实施例。
[0107]
实施例1
[0108]
下面以一个变构蛋白质(pdb id：2xem)为例，来介绍allosests方法预测蛋白质变构口袋的实施过程。已知该蛋白质的三维结构，通过实施该方法可获得其变构口袋的位置及排名。预测过程需在linux系统下完成。
[0109]
(1)准备工作
[0110]
在运行蛋白质变构口袋预测程序前，需确保已安装fpocket、python 3.6、numpy、pandas、sklearn、joblib和prody工具包。蛋白质口袋内氨基酸残基的二级结构类型可使用dssp工具完成。此外，还需使用psi-blast程序(https://blast.ncbi.nlm.nih.gov/blast.cgi)计算蛋白质序列位置特异性打分矩阵pssm。将处理好的蛋白质三维结构及预先计算得到的特征数据存放在allosests根目录下。
[0111]
(2)预测蛋白质变构口袋
[0112]
在完成上述准备工作后，可运行变构口袋预测程序predict.py。首先将处理后的2xem.pdb文件、predict.py和utils.py程序、标准pssm和dssp文件放在同一路径下。
[0113]
执行以下命令，运行allosests程序来预测蛋白质变构口袋：
[0114]
…
]$python predict.py
‑‑
pdbid[pdbid]
‑‑
chain[chain]
[0115]
例：
[0116]
…
]$python predict.py
‑‑
pdbid 2xem
‑‑
chain b
[0117]
其中，[pdbid]和[chain]为输入参数，[pdbid]为蛋白质pdb文件名称，[chain]为所要预测的功能链。
[0118]
执行上述命令后，allosests程序会查找蛋白质表面口袋、提取所需要的口袋特征，最终给出口袋的位置和预测概率排名(打印在bash命令界面，如附图4所示)。
[0119]
实施例2
[0120]
将本发明构建的allosests模型在两个独立测试集上进行测试。对测试集中的每一个变构蛋白质的计算过程同实施例1，这里只对结果作说明。将allosests与目前发展的三种蛋白质变构口袋预测方法pars、allosite和allositepro在测试集1上进行预测性能比较，其结果见表6；allosests与passer2.0在测试集1和2上进行预测排名比较，结果见表7。
[0121]
表6 allosests与其他方法在测试集1上的预测性能比较
[0122][0123]
从四种方法在测试集1上的预测性能来看(表6)，本发明提出的算法allosests具有更好的预测蛋白质变构口袋的能力，整体预测性能显著优于另三种方法，其敏感性sen为0.750，特异性spe为0.932，精度pre为0.474，和马修斯相关系数mcc为0.556。
[0124]
表7 allosests与其他方法在测试集1和2上预测口袋排名的比较
[0125][0126]
top 1,top 2,top 3:排名前1、2、3位。
[0127]
从两种方法在测试集1和2上的预测排名来看(表7)，虽然allosests将较少的蛋白质变构位点(测试集1中的45.8％和测试集2中的60.7％)排在top 1，但在top 2和top 3上实现了更高的百分比，即在测试集1上将75.0％和83.3％的蛋白质的真实变构口袋排在了top 2和top 3；在测试集2上将78.6％和89.3％的蛋白质的真实变构口袋排在了top 2和top 3。
[0128]
以上结果说明，allosests方法中新增加的口袋特征，包括二级结构信息、动力学转移熵和空间协同进化信息，以及实施的基于多份训练子集的集成模型构建策略确实能够改进模型对蛋白质变构口袋的预测性能，它们的引入对蛋白质变构口袋预测有较好的贡献。