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

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

背景技术：

1、变构是蛋白质发挥功能的重要调控手段，通常是指调节分子(变构剂)结合于蛋白质活性部位(正构位点)以外的其它部位(变构位点)，引起活性部位结构或动力学性质的改变，从而影响蛋白质功能的现象[1-3]。比较正构与变构位点的空间和拓扑结构显示变构位点的结构更具多样性，这使得变构调节剂具有更高选择性和更少毒副作用，这些优势使得变构剂设计在药物研发中具有更重要的意义[4-6]。变构位点的发现是变构药物设计的前提，如何准确高效地识别蛋白质变构位点是结构生物学家和药物设计者普遍关心的问题。

2、目前，实验上获得变构分子/变构位点多是由高通量筛选实验偶然所得[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]算法来解决特征工程、模型选择和超参调优等耗时费力的问题。

3、最新研究表明，蛋白质变构活性在本质上是熵的效应，它依赖于由残基间协同涨落引起的信息从变构位点向正构位点的传递[25-26]。近来，由schreiber[27]提出的转移熵概念，由于其与信息流和变构交流有关，得到了广泛关注。转移熵中引入了残基间时间延迟运动相关性信息，因此它可以呈现并量化残基间驱动与被驱动的关系，从而揭示哪些位点可以被操控来影响蛋白质正构位点的活性，这些位点很可能成为蛋白质的变构位点。另外，蛋白质功能位点残基，如配体结合和催化位点，往往不是独立发挥作用，而是彼此协同发挥功能[28-30]，蛋白质变构口袋残基也一样。这种协同性必然会体现在残基进化上，从而呈现出一定的空间协同进化特点。最近，我们提出了一种新的进化信息编码方法，即基于空间近邻的位置特异性打分矩阵(spacial neighbor-based position-specific scoring matrix,snb-pssm)，它在编码残基进化信息时通过近邻平均的方式在一定程度上考虑了残基空间协同进化的特点[31]，在预测蛋白质上结合rna的位点中取得了好的效果[32]。以上是从特征方面提出了可能的改进策略，从模型的角度看，集成模型由于综合了多个子模型的预测结果往往比单一模型具有更好的泛化能力。在此，我们提出了基于多份训练子集的集成模型构建策略。它将大类(包含样本数相对较多的阳性或阴性类)样本拆成数份、并将每一份与小类(包含样本数相对较少的阳性或阴性类)样本组合、产生多份训练子集、选择适合多份训练子集的最优特征组合、之后构建多个子模型，并将子模型的预测结果进行整合得到一个集成模型。该方法不仅可以增强模型的泛化能力而且在一定程度上解决了样本不平衡的问题。基于以上分析，我们在蛋白质变构位点预测中考虑了转移熵和残基协同进化信息，并采用基于多份训练子集的集成模型构建策略，提出了考虑转移熵和空间近邻进化信息的基于集成模型的蛋白质变构位点预测方法。

4、参考文献

5、[1]greener j g,sternberg m j.structure-based prediction of proteinallostery[j].curr opin struct biol.2018,50:1-8.

6、[2]liu j,nussinov r.allostery:an overview of its history,concepts,methods,and applications[j].plos comput biol.2016,12(6):e1004966.

7、[3]zha j,li m,kong r,et al.explaining and predicting allostery withallosteric database and modern analytical techniques[j].j mol biol.2022,434(17):167481.

8、[4]lu s,he x,ni d,et al.allosteric modulator discovery:fromserendipity to structure-based design[j].j med chem.2019,62(14):6405-6421.

9、[5]guarnera e,berezovsky i n.allosteric sites:remote control inregulation of protein activity[j].curr opin struct biol.2016,37:1-8.

10、[6]cheng x,jiang h.allostery in drug development[j].adv exp medbiol.2019,1163:1-23.

11、[7]jiang y.,kalodimos c.g.nmr studies of largeproteins.j.mol.biol.2017；429:2667–2676.

12、[8]xiao s,verkhivker g m,tao p.machine learning and protein allostery[j].trends in biochemical sciences.2022.

13、[9]gulati s,palczewski k,engel a,et al.cryo-em structure ofphosphodiesterase 6 reveals insights into the allosteric regulation of type iphosphodiesterases[j].sci adv.2019,5(2):v4322.

14、[10]qi y,wang q,tang b,et al.identifying allosteric binding sites inproteins with a two-statemodel for novel allosteric effector discovery[j].jchem theory comput.2012,8(8):2962-2971.

15、[11]weinkam p,pons j,sali a.structure-based model of allosterypredicts coupling between distant sites[j].proc natl acad sci u s a.2012,109(13):4875-4880.

16、[12]goncearenco a,mitternacht s,yong t,et al.spacer:server forpredicting allosteric communication and effects of regulation[j].nucleicacids res.2013,41(web server issue):w266-w272.

17、[13]panjkovich a,daura x.pars:a web server for the prediction ofprotein allosteric and regulatory sites[j].bioinformatics.2014,30(9):1314-1315.

18、[14]ma x,meng h,lai l.motions of allosteric and orthosteric ligand-binding sites in proteins are highly correlated[j].j chem inf model.2016,56(9):1725-1733.

19、[15]suel g m,lockless s w,wall m a,et al.evolutionarily conservednetworks of residues mediate allosteric communication in proteins[j].natstruct biol.2003,10(1):59-69.

20、[16]reynolds k a,mclaughlin r n,ranganathan r.hot spots forallosteric regulation on protein surfaces[j].cell.2011,147(7):1564-1575.

21、[17]amor b r,schaub m t,yaliraki s n,et al.prediction of allostericsites and mediating interactions through bond-to-bond propensities[j].natcommun.2016,7:12477.

22、[18]wang j,jain a,mcdonald l r,et al.mapping allostericcommunications within individual proteins[j].nat commun.2020,11(1):3862.

23、[19]huang w,lu s,huang z,et al.allosite:a method for predictingallosteric sites[j].bioinformatics.2013,29(18):2357-2359.

24、[20]le guilloux v,schmidtke p,tuffery p.fpocket:an open sourceplatform for ligand pocket detection,bmc bioimformatics 2009；10:168.

25、[21]song k,liu x,huang w,et al.improved method for the identificationand validation of allosteric sites[j].j chem inf model.2017,57(9):2358-2363.

26、[22]xiao s,tian h,tao p.passer2.0:accurate prediction of proteinallosteric sites through automated machine learning[j].front mol biosci.2022,9:879251.

27、[23]jin h.,song q.,x h.auto-keras:an efficient neural architecturesearch system:in proceedings of the 25th acm sigkdd international conferenceon knowledge discovery&data mining[z].20191946-1956.

28、[24]erickson n,mueller j,shirkov a,et al.autogluon-tabular:robust andaccurate automl for structured data[j].arxiv prepr.2020,arxiv:2003.06505.

29、[25]tsai cj,nussinov r.a unified view of"how allostery works",ploscomputational biology 2014；10:e1003394.

30、[26]motlagh h n,wrabl j o,li j,et al.the ensemble nature of allostery[j].nature.2014,508(7496):331-339.

31、[27]schreiber t.measuring information transfer[j].phys rev lett.2000,85(2):461-464.

32、[28]yang z,deng x,liu y,et al.analyses on clustering of the conservedresidues at protein-rnainterfaces and its application in binding siteidentification[j].bmc bioinformatics.2020,21(1):57.

33、

34、[29]guharoy m,chakrabarti p.conserved residue clusters at protein-protein interfaces and their usein binding site identification[j].bmcbioinformatics.2010,11:286.

35、[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.

36、[31]liu y,gong w,yang z,et al.snb-pssm:a spatial neighbor-based pssmused forprotein-rna binding site prediction[j].j mol recognit.2021,34(6):e2887.

37、[32]liu y,gong w,zhao y,et al.aprbind:protein-rna interfaceprediction by combiningsequence and i-tasser model-based structural featureslearned with convolutional neuralnetworks[j].bioinformatics.2021,37(7):937-942.

技术实现思路

1、本发明的目的是在已知变构蛋白质结构的情况下，预测其变构位点，为蛋白质变构位点识别和变构药物设计提供重要信息。方法的特点在于：①查找蛋白质表面潜在的变构口袋，并提取口袋特征。除口袋的理化性质外，本发明还提取了转移熵和残基协同进化信息。结果显示这些新加入的特征对变构位点预测有好的贡献，在特征选择中被选入最优特征组合。②本发明提出了基于多份训练子集的集成模型构建策略，即拆分大类样本、产生多份训练子集、选择适合多份训练子集的最优特征组合、构建基于多个子模型的集成模型来进行变构位点预测。该策略由于利用了训练集中所有的大类样本，因此它除了有利于模型的泛化能力提高外，还附带解决了训练集中变构位点正负样本极度不平衡的问题。③本发明提出了基于多份训练子集的集成模型构建策略，并首次将转移熵和空间近邻进化信息用于蛋白质变构位点预测，其中空间近邻进化信息是我们之前开发的(jmol recognit.2021,34(6):e2887)，可以很好地考虑氨基酸残基空间近邻的协同进化性。

2、本发明考虑转移熵和空间近邻进化信息的基于集成模型的蛋白质变构位点预测方法，是一种基于转移熵和空间近邻进化信息的集成分类模型来对蛋白质变构位点进行预测的方法(ensemble-based machine learning approach with transfer entropy andsnb-pssm considered for protein allosteric site prediction,allosests)，包括四个步骤：一是查找蛋白质表面潜在的变构口袋，二是提取口袋特征，三是产生多份训练子集并筛选获得最优特征组合，四是构建基于多个子模型的集成模型以预测蛋白质变构位点(如附图1)。

3、建立模型前，首先构建了训练集和测试集：

4、收集来自蛋白质数据库(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 cheminf 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、链号、残基数和口袋样本数；

5、表1训练集中的127个变构蛋白质及其口袋样本数

6、

7、     

8、     

9、     

10、

11、a/b:变构口袋数/口袋总数。

12、表2测试集1中的24个变构蛋白质及其口袋样本数

13、

14、

15、a/b:变构口袋数/口袋总数

16、表3测试集2中的28个变构蛋白质及其口袋样本数

17、

18、

19、a/b:变构口袋数/口袋总数

20、下面对蛋白质变构位点预测方法的四个步骤进行描述。

21、步骤1：查找蛋白质表面潜在的变构口袋

22、针对选定的变构蛋白质，本发明保留带有变构残基的功能链，去除水分子和配体，将其用于查找蛋白质表面潜在的变构口袋，这些变构口袋被作为样本数据；样本数据中与变构调节剂原子接触最多的口袋被标记为阳性样本，其余被标记为阴性样本；对于含有127个蛋白质的训练数据集，查找到2245个口袋，其中有147个阳性样本和2098个阴性样本；对于两个测试集，即包含24个蛋白质的测试集1(test set 1)和包含28个蛋白质的测试集2(test set 2)，分别查找到24和28个阳性样本，295和347个阴性样本；

23、这里，使用高效开源的fpocket(bmc bioimformatics 2009,10:168)工具查找蛋白质表面潜在的变构口袋，采用默认的参数设置；对于查找到的蛋白质表面口袋，fpocket提供了19个物理化学特征，如表4所列。这些口袋如前所述被分成训练集和测试集用于预测模型的训练和测试；

24、表4fpocket计算的口袋理化性质特征

25、

26、

27、步骤2：提取蛋白质口袋特征

28、对查找到的蛋白质表面潜在变构口袋，除19个理化性质特征(表4)外，还提取了16个表征口袋二级结构信息的特征、20个表征空间协同进化信息的特征，和1个表征动力学转移熵信息的特征。每一个蛋白质表面口袋对应的后三类特征的提取具体描述如下：

29、口袋二级结构信息的特征：是根据其中残基的8种二级结构类型来提取的，包括口袋中氨基酸残基在每种二级结构类型中出现的数量，和该数量占口袋内全部残基数量的百分比。8种二级结构类型包括：α-helix(h)、β-bridge(b)、β-ladder(e)、310-helix(g)、π-helix(i)、turn(t)、bend(s)和不确定类型(m)，二级结构类型使用dssp工具(biopolymers.1983,22(12):2577-2637)计算；

30、空间协同进化信息的特征：我们之前(j mol recognit.2021,34(6):e2887)提出的基于空间近邻的位置特异性打分矩阵(snb-pssm)被用于提取蛋白质口袋的进化信息。首先，对目标蛋白质序列，用psi-blast工具(nucleic acids res.1997,25(17):3389-3402)在非冗余蛋白质序列数据库中循环查找其同源序列，以e-value＝0.001为阈值迭代运行三次，生成标准位置特异性打分矩阵(position-specific scoring matrix,pssm)。pssm给出了同源序列比对的结果，其中包含进化信息；对于含有n个残基的蛋白序列，pssm矩阵的大小为20×n，每个残基位置的进化信息被编码在20维的向量中；然后，为了考虑目标残基空间近邻的进化信息，对其原始进化信息进行基于空间近邻的平滑处理，使其进化信息为：与该残基ca原子距离以内的所有残基进化信息的平均值(见附图2的一个例子)。最后，口袋内所有残基进化信息的平均值被作为描述口袋的空间协同进化信息特征，该特征为一个20维的向量。

31、动力学转移熵信息的特征：动态高斯网络模型(dynamic gaussian networkmodel,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+τ)值不确定性的减少量：

32、ti→j(τ)＝s(△rj(t+τ)|△rj(t))-s(△rj(t+τ)|△ri(t),△rj(t)) (1)

33、这里，s(δrj(t+τ)|δrj(t))是给定δrj(t)时，δrj(t+τ)的条件熵；s(δrj(t+τ)|δri(t),δrj(t))是给定δri(t)和δrj(t)时，δrj(t+τ)的条件熵；根据香农熵，等式(1)中的两项为：

34、s(△rj(t+τ)|△rj(t))＝-<ln p(△rj(0),△rj(τ))>+<ln p(△rj(0))> (2)

35、

36、将等式(2)和(3)代入(1)得到：

37、

38、其中，<>表示系综平均，p()表示概率。

39、在gnm模型下(以cα原子为节点，截断半径为时间延迟交叉相关和转移熵分别为：

40、

41、

42、这里，和分别是第k个运动模式下残基i与j的交叉相关性、残基j的自相关和残基i的自相关，τ0为特征时间(characteristic time)，这里采用τ/τ0＝5的时间延迟。运动模式由gnm下蛋白质体系的基尔霍夫矩阵获得，基尔霍夫矩阵为：

43、

44、其中，为截断半径，rij表示两个节点间的距离。基尔霍夫矩阵的伪逆矩阵为：

45、

46、其中，k为运动模式数，λk为对应第k个运动模式的特征值，μk为对应第k个运动模式的特征向量。

47、从残基i到j的净熵转移(net entropy transfer,net)是从i到j的熵转移与从j到i的熵转移之差：

48、neti→j＝ti→j(τ)-tj→i(τ) (9)

49、从残基i到蛋白质其他残基的净熵转移neti→e是上式对残基j求和。为方便比较，对neti→e进行归一化处理。口袋的净熵转移被定义为口袋中所有残基净熵转移的平均值，该值若大于0，则口袋倾向于发出信号，否则倾向于接收信号；口袋的净熵转移即为口袋动力学转移熵信息的特征。

50、步骤3：产生多份训练子集并筛选获得最优特征组合

51、在变构位点预测中，步骤1中对应的正负样本数量是极不平衡的，可达到1:10。如此高的不平衡将严重影响分类器模型的预测性能，造成对大类样本的过拟合。通常，过采样或欠采样技术被用于处理不平衡的数据。过采样是通过随机复制或生成小类样本，而欠采样则是通过删除部分大类样本，来达到两类数据的平衡。两种技术分别可能会产生不合理的样本或造成一定程度信息的损失。为此，我们提出了一种基于多份训练子集的集成模型构建策略来克服样本的不平衡问题。首先，将训练集中的大类样本随机分为10个子集(子集中样本数与小类样本数相当)，并将它们分别与小类样本组合成10份训练数据集；接下来，对每一份训练数据集，使用前向特征选择算法(sequential forward selection，sfs)从预先提取的特征(步骤2中的56种特征)中筛选出最优特征子集；这里，我们通过对每一份训练数据集进行5次10-fold交叉验证来获得最优特征子集，然后将10个最优特征子集取并集获得一个特征子集；对该特征子集，使用最大相关最小冗余算法(maximum correlationminimum redundancy,mrmr)进一步筛选获得最终的最优特征组合；mrmr能够保证所选特征与输出变量间的最大相关，同时保证所选特征间的最小冗余；表5列出了最终筛选出的最优特征组合；

52、表5特征选择筛选的最优特征组合

53、

54、

55、步骤4：构建基于多个子模型的集成模型以预测蛋白质变构位点

56、针对每一份训练数据集，建立相应的支持向量机(support vector machine,svm)分类子模型，通过5-fold交叉验证来训练模型；然后对10个子模型的预测结果(即样本是变构口袋的概率)进行平均整合，构建集成分类模型，来预测口袋是否为蛋白质的变构口袋(如附图3的一个例子)。输入数据为表5列出的最优特征组合数据，输出为预测的蛋白质变构口袋的位置和概率排名；

57、allosests模型在10份训练集上均使用5-fold交叉验证训练和调参(通过最大化马修斯相关系数mcc)，并在独立测试集1和2上进行测试；模型的评价指标包括敏感性(sensitivity,sen)、特异性(specificity,spe)，精度(precision,pre)和马修斯相关系数(matthews correlation coefficient,mcc)，其公式如下：

58、

59、

60、

61、

62、其中，tp(true positive)、fp(false positive)，tn(true negative)和fn(falsenegative)分别为真阳性、假阳性，真阴性和假阴性样本数量，它们通过比较每个口袋的预测情况与真实情况来确定。对于模型的测试，我们将测试蛋白质直接输入到集成模型。

63、本发明提出了一种考虑转移熵和空间近邻进化信息的基于集成模型的蛋白质变构位点预测方法allosests。该方法能够比较快速地预测蛋白质变构口袋的位置并进行排名，在预测性能上比现有的其他预测方法具有更好的评价指标。与其他预测方法相比，本发明不仅考虑了蛋白质口袋的理化性质和结构特征，还考虑了动力学转移熵信息和残基协同进化信息，是一个基于多个子模型的集成分类器。allosests在蛋白质变构口袋预测中有较好性能，可用于蛋白质变构药物靶点识别和变构药物设计的研究。