CN111951884A - 蛋白质小分子结合口袋上关键柔性氨基酸的识别方法 - Google Patents

Abstract

本发明涉及一种蛋白质小分子结合口袋上关键柔性氨基酸的识别方法，包括以下步骤：步骤1、在蛋白质PDB结构中，通过95％的序列相似性将蛋白质的结构分类为不同的亚型簇；步骤2、明确的生物学配体，将在同一结构亚型簇中所有蛋白结构重叠，并准确鉴定蛋白上小分子结合口袋；步骤3、将已知3D结构的同一结构亚型簇中的相应结合位点，根据小分子结合位置划分为位点结合亚型；步骤4、通过计算理化和结构参数差异表征因子来研究蛋白质柔性的性质；步骤5、用形成蛋白质结合口袋的系统的最高评分来识别关键的柔性残基。设计了一套统一性、无偏见的研究流程体系，开发了一套系统性的数学统计公式来确定不同结构模型相同位置上的结合位点氨基酸的结构柔性。

Description

蛋白质小分子结合口袋上关键柔性氨基酸的识别方法

技术领域

本发明涉及计算生物技术领域，尤其涉及一种蛋白质小分子结合口袋上关键柔性氨基酸的识别方法。

背景技术

蛋白质是组成生命的物质基础，蛋白三维结构是体现自身生物化学功能和细胞功能的基本条件。蛋白质在执行生物功能过程中，不断改变自身构象并与周围环境产生理化和生物反应。这种结构柔性特征，与蛋白自身结构功能域、分子识别与结合机制和催化活性等生理特性息息相关，对于蛋白质正常生理功能和生物的发挥有着极其重要的作用。

目前，基于分子对接原理的药物虚拟筛选技术，成为药物研发的重要研究工具。蛋白质上组成小分子结合口袋的氨基酸，其结构可塑性(plasticity)特别是部分柔性氨基酸残基，对分子对接结果精度有着决定性影响。因此，鉴定小分子结合口袋中关键柔性氨基酸位点，对于基于蛋白可塑性特征的分子对接技术在药物研发等方面，具有极为重要的指导意义。

随着结构生物学和计算生物学的飞速发展，RCSB PDB数据库已解析出>16万蛋白晶体结构，极大地促进了蛋白质结构柔性状态的结构研究。目前，已经开发了几种基于现有晶体结构数据预测蛋白构象变换的经验性方法。MolMovDB数据库涉及蛋白质和RNA上从小环、结构域到整个亚基二级结构柔性预测。DynDom程序的核心是根据动态域、域间螺旋轴和域间弯曲区域构象变化的结果而开发的。PCDB数据库通过对非冗余性、各种类蛋白质上的结构域间最大RMSD值进行计算和比较分析，获取各二级结构域上柔性状况。根据分子动力学方法模拟轨迹，计算每个氨基酸的多个RMSD数据，是CABS-flex服务器探测研究蛋白质柔性状态的主要方法。CoDNaS首先收集同种蛋白质非冗余性3D结构集合，并根据C-alphaRMSD作为距离度量来探索构象异构体之间的结构差异。

除上述数据库外，针对的特定蛋白质亚型的研究也有很多。但是，系统性探索蛋白整体或局部区域构象柔性特征，并在氨基酸水平上确定结合位点和在原子水平预测关键柔性残基性质的研究并不多。

发明内容

(一)要解决的技术问题

鉴于现有技术的上述缺点、不足，本发明提供一种蛋白质小分子结合口袋上关键柔性氨基酸的识别方法，来确定不同结构模型相同位置上的结合位点氨基酸的结构柔性。

(二)技术方案

为了达到上述目的，本发明采用的主要技术方案包括：

本发明一种蛋白质小分子结合口袋上关键柔性氨基酸的识别方法，包括以下步骤：

步骤1、在蛋白质PDB结构中，通过95％的序列相似性将蛋白质的结构分类为不同的亚型簇；

步骤2、明确的生物学配体，将在同一结构亚型簇中所有蛋白结构重叠，并准确鉴定蛋白上小分子结合口袋；

步骤3、将已知3D结构的同一结构亚型簇中的相应结合位点，根据小分子结合位置划分为位点结合亚型；

步骤4、通过计算理化和结构参数差异表征因子来研究蛋白质柔性的性质；

步骤5、用形成蛋白质结合口袋的系统的最高评分来识别关键的柔性残基。

可选地，所述步骤1中，采用>95％序列相似性阈值，利用Cd-hit软件对共计133976个PDB ID蛋白结构数据集，剔除冗余序列，从而得到总共41754个亚型结构簇和410079个非冗余链序列；

其中多蛋白链簇：具备多条蛋白链(至少2条链)的集合簇，删除不符合条件的6591个数据簇集，最终得到35963个结构数据集。

可选地，生物配体与蛋白质上邻近氨基酸残基的接触率，被视为配体与周围环境相互作用程度的几何量度；如果两个重原子的中心位于不超过范德华半径之和＞20％的距离处，则认为它们之间存在相互接触；

溶剂可及化区域(ASA)指小分子配体暴露于周围水溶液环境中的程度；DMS软件被用来计算41754个蛋白质簇中、不同结构状态下的生物活性化学分子的溶剂可及化区表面积值，公式如下：

ASA ratio＝1-ASA(complex)/ASA(apo)

如果小分子配体能够同时满足接触率≥60％和ASA率≤40％的阈值标准，则表明该分子与周围氨基酸具有良好的结合状态；在28853蛋白质结构簇中，当考虑至少四个氨基酸残基与小分子接触时，小分子配体并没有位于合理的位置；此外，当研究单个簇成员时，在所有6310个簇中对98693链的53580个PDB ID蛋白上小分子结合口袋位于合理位置。

可选地，利用ClustW程序，对同一聚类中不同蛋白链进行多个序列比对；基于序列比对信息，将蛋白中位于同一位置、氨基酸不同的非保守性位点设定为单个序列中的突变位点；

利用DSSP程序定义了二级结构和几何特征的对蛋白质数据格式，用于校准5996个蛋白质簇和95007链的二级结构和氨基酸位点，以获得氨基酸序列和位置编码信息；

氨基酸序列归一化：根据序列相似性保守程度和位点校准信息，对蛋白簇中所有链的氨基酸序列进行归一化和编码，以使所有蛋白质在相同位置具有相同的氨基酸编码；

突变：结合氨基酸序列归一化编码信息的结构和几何特征，可以获得每个蛋白质链相同的二级结构区域信息和氨基酸序列突变位点信息。

可选地，基于几何构型，通过搜索与配体相邻氨基酸的方法进行小分子结合口袋的预测；为最大化囊括所有氨基酸，将每条链中距离生物活性配体分子重原子

以内的所有残基设定为分子结合袋；

来自同一蛋白簇、不同的结构链中，形成结合袋的残基彼此之间存在显着差异；与蛋白4GTZ链A相比，由于小分子配体分子不同，4GTY蛋白A链中存在5个新的形成小分子配体的氨基酸，以圆圈标记；此外，如果同一链中不同配体上重原子的距离小于

则这些配体被视为一个整体的化合分子；而且，若同一蛋白质链上存在不同的结合口袋区域，并且口袋之间不同区域的重叠率高于>60％，这些不同结合位点被整合为同一个口袋。

可选地，对于各亚型结构簇，模板结构选时，优先级高的元素选择设定为单个模板结构：

拥有生物小分子配体>无突变>氨基酸序列长度最长>分辨率最高。

可选地，使用具有相似拓扑结构的已知模板，在蛋白结构表面上探索口袋或空腔；基于模板的方法，通过全局拓扑搜索同源结构；而后，假定的结合位点可以在蛋白结构叠合后进行坐标转化；使用软件UCSF Chimera根据结合口袋主链重原子叠加同一簇中的所有蛋白质链，以获得蛋白质的叠加簇；每条链、相同位置处的结合位点氨基酸被收集为一个集合体，其被称为配体位点整合体。

可选地，RMSD:重原子位置均方根偏差(RMSD)值，是测量不同蛋白质体系上原子间平均距离的一种度量形式，在此计算了模板结构和和目标结构之间每个结合位点氨基酸上所有重原子、侧链和骨架主链RMSD值；包括ASP，GLU，PHE和TYR氨基酸在内的侧链末端对称原子被视为同一原子；

C-factor：在典型的蛋白质晶体分辨率参数中，高温度因子反映了重原子的低经验电子密度指数，反之亦然；温度因子值，被存储在各PDB格式文件内ATOM和HETATM坐标信息中，因此B因子可以反应氨基酸原子的真实静态或动态迁移率；但是，由于结晶过程中，应用了不同的结构优化程序，在不同已知结构中测量B因子的规模可能不同；因此，在本研究中，使用标准化数据来代替原始数据，以比较不同蛋白质结构链温度B-factor因子之间差异；原子归一化的B因子值计算如下：

Average B-factor value

Sample standard deviation of B-factor value：

Each atom normahzed B-factor value：

B_{normalized-single}＝(B-factor-B_average)/B_{normalized-all}

不同结构与模板结构之间K277氨基酸上的B因子和重原子均方根偏差RMSD差异较大为其侧链结构，骨架重原子差异较小；氨基酸P305上所有重原子、K277的4骨架主链原子中仅观察到微小的结构波动；相较于模板结构，其它蛋白链上K277氨基酸上的侧链原子间的温度因子B-factor和Rmsd差异较为明显的，表明此处氨基酸较为柔性。

可选地，溶剂可及化表面比率，是通过计算不同结构状态下、相同氨基酸位点的溶剂可及化表面积比值，利用DMS软件计算在每条链中形成配体结合位点的所有氨基酸的ASA值。

可选地，计算和比较空间位阻：对各配体分子与所有整体口袋中所有氨基酸之间的空间位阻进行计算和比较分析；设定各重原子之间接触参数，并利用UCSF Chimera软件计算各原子间空间位阻；若两重原子之间范德华半径之和减去它们之间的距离并减去可能存在的氢键对的余量，两重原子之间的结构被视为结构重叠或空间位阻：

Overlap_ij＝R_VDWi+R_VDWj-d_ij-allowance_ij

Overlap Cutoff＝1.H-bond Allowance＝0.6

为了检测配体分子与氨基酸之间空间位阻，采用阈值标准临界值为

浮动值为

(默认空间位阻参数标准分别为0.6和

)；较大的正截距值，表征空间位阻较为显著；而负截距值，表明两重原子之间并不存在空间位阻；不同结构上、相同位置的氨基酸之间若存在较大空间位阻，表明此处氨基酸构象变化大；随着残基K277的构象变化，蛋白质A链上的氨基酸K277与蛋白质C上的配体存在明显的空间位阻，表明K227氨基酸位点区域柔性变化显著；

鉴定突变体的突变位点：

突变引起的轻微结构扰动，会导致整个蛋白质发生变化；在同一蛋白簇中，通过多序列比对和基于序列比对获取的位点信息，鉴定各蛋白质链的突变位点信息；蛋白质序列中突变频率高，表明蛋白质在此处柔性较大；

根据突变前后氨基酸的几何片段和理化性质的变化，将突变位点分为不同等级，并赋予不同的权重系数值；其中，突变前后氨基酸若具有完全相反的物理和化学特性，将对该位点赋予较高的权重系数，为突变之前和之后具有完全相反的物理和化学特性的氨基酸分配更高的加权因子。

可选地，计算二面角：每个残基的二面角主要集中在侧链重原子上，并通过UCSFChimera软件进行计算；

确定关键柔性氨基酸位点

在每一个蛋白结构聚类中，赋予可表征蛋白质氨基酸柔性的权重因子不同系数，公式如下：

可选地，突变信息和蛋白-小分子间空间位阻，被认为氨基酸柔性的最主要位点，都被赋予1.0的重要系数；

反应不同氨基酸位相间差异的二面角、RMSD值、溶剂可及化表面积值，被赋予0.9的权重系数，B-factor值权重系数为0.8；

根据不同权重系数表征因子计算公式，最后得到综合性的打分，并对组成配体结合位点的所有氨基酸进行打分和排名，排名靠前的氨基酸残基表明该位点结构位相差异较大。

(三)有益效果

基于RCSB PDB数据库中共计156800个结构，设计了一套统一性、无偏见的研究流程体系，并通过计算6个综合“表征因子”(RMSD值，B因子，突变位点，溶剂可及化表面积，空间位阻和二面角)以标记蛋白小分子结合位点上氨基酸残基的结构变化信息。最后，本发明开发了一套系统性的数学统计公式来确定不同结构模型相同位置上的结合位点氨基酸的结构柔性。

本发明不仅在蛋白质整体结构水平、以二级结构为主的局部区域水平上，系统性的探索了蛋白整体或构象柔性特征，并且在氨基酸水平上确定了小分子结合位点，并基于原子水平对氨基酸侧链和骨架原子进行了系统性分析，得到了信息更加丰富的关键氨基酸位点柔性特征数据。该方法所鉴定的关键柔性氨基酸特征数据，在分子对接技术中，能够有效减少氨基酸刚性位相空间阻碍，提高分子对接精度和药物分子虚拟筛选精确率。本发明对于基于蛋白可塑性特征的分子对接技术在药物研发等方面，具有极为重要的意义。

在真实的数据集上的实验结果表明，与多个现有方法相比，本方法具有更高的准确度。

附图说明

图1为蛋白质小分子结合口袋上关键柔性氨基酸鉴别流程图；

图2为不同PDB ID编码晶体结构4GTZ和4GTY上A链小分子结合口袋图A图；

图3为不同PDB ID编码晶体结构4GTZ和4GTY上A链小分子结合口袋图B图；

图4为蛋白质内部特异性亚型结构簇分类二维示意图；

图5为组成亚型聚类结构簇内配体结合口袋“整合体”周围氨基酸残基；

图6为配体结合口袋周围氨基酸结构差异化分析A示意图；

图7为配体结合口袋周围氨基酸结构差异化分析B示意图；

图8为配体结合口袋周围氨基酸结构差异化分析C示意图；

图9为配体结合口袋周围氨基酸结构差异化分析D示意图；

图10为配体结合口袋周围氨基酸结构差异化分析E示意图；

图11为配体结合口袋周围氨基酸结构差异化分析F示意图；

图12为配体结合口袋周围氨基酸结构差异化分析G示意图；

图13为各氨基酸突变位点被赋予不同权重系数值二维示意图；

图14为突变赋予权重系数图。

具体实施方式

为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。

为了更好的理解上述技术方案，下面将参照附图更详细地描述本发明的示例性实施例。虽然附图中显示了本发明的示例性实施例，然而应当理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更清楚、透彻地理解本发明，并且能够将本发明的范围完整的传达给本领域的技术人员。

具体实施例描述部分：

基于上述原理，本发明将从蛋白质PDB结构中(1)通过95％的序列相似性将蛋白质的结构分类为不同的结构亚型簇；(2)明确的生物学配体，将在同一结构亚型簇中所有蛋白结构重叠，并准确鉴定蛋白上小分子结合口袋；(3)将已知3D结构的同一结构亚型簇中的相应结合位点，根据小分子结合位置划分为位点结合亚型；(4)通过计算理化和结构参数差异“表征因子”来研究蛋白质柔性的性质；(5)用形成蛋白质结合口袋的系统的最高评分来识别关键的柔性残基。通过该分析，鉴定出了6993个蛋白亚型簇，即小分子配体子口袋集合簇，共计79244个PDB ID蛋白和135760条结构链。

材料和方法

数据集

蛋白质：蛋白质晶体数据库(截止日期：2019-10-23)中，共计含有总数为156800个PDB ID蛋白质结构、148303个多肽链生物单位结构。所有数据均从ftp://ftp.wwpdb.org网站下载。其中，首先将核磁共振NMR法获取的晶体结构排除。其次，仅仅包含RNA、DNA和非标准残基类晶体结构也被排除在外。与PDB ID蛋白结构相比，具有更多单元链的多肽链结构优先考虑。蛋白质氨基酸序列：每种蛋白质的完整氨基酸序列，被临时存储以用于后续的聚类分析。肽：截断少于30个氨基酸的独立链结构被定义为肽。

生物配体：提取每个结构中带有HETATM注释的化学分子。但是，RCSB PDB数据库的化学结构通常包含溶剂类物质如水分子和离子(Na⁺，Cl^-和Zn²⁺等)需要排除在外；尤其是部分结晶条件的非生物活性分子(缓冲液和辅助因子)也要排除在外。此外，与特异靶标具有较好结合能力的共价键类分子配体也不考虑在内。因此，将蛋白结构中含有的具备生物学条件的化学分子、拮抗剂、抑制剂和药物分子，归类为生物活性配体。参照其它研究文献，在本申请方案中，生物活性配体参照标准如下：(a)非氢重原子总数：>5；(b)碳原子数：>3；(c)重原子的分子量：60～100道尔顿。

序列聚类

聚类：Cd-hit是目前最为广泛的、用于序列聚类以剔除冗余蛋白序列的程序之一。采用>95％序列相似性阈值，利用Cd-hit软件对共计133976个PDB ID蛋白结构数据集，剔除冗余序列，从而得到总共41754个亚型结构簇和410079个非冗余链序列。

本申请方案旨在鉴定不同蛋白链之间、相同结合口袋位置上重要的氨基酸残基的构象变化。

多蛋白链簇：具备多条蛋白链(至少2条链)的集合簇，删除不符合条件的6591个数据簇集，最终得到35963个结构数据集。非空簇：小分子化合物是定位蛋白质上活性结合口袋的最重要指标。因此，选择具有至少1种生物配体的非空蛋白链的最终数据集，用于结合位点分析。

原子接触率和溶剂可及表面积

生物配体与蛋白质上邻近氨基酸残基的接触率，被视为配体与周围环境相互作用程度的几何量度。如果两个重原子的中心位于不超过范德华半径之和>20％的距离处，则认为它们之间存在相互接触。

溶剂可及化区域(ASA)指小分子配体暴露于周围水溶液环境中的程度。DMS软件被用来计算41754个蛋白质簇中、不同结构状态下的生物活性化学分子的溶剂可及化区表面积值，公式如下：

ASA ratio＝1-ASA(complex)/ASA(apo)

如果小分子配体能够同时满足接触率≥60％和ASA率≤40％的阈值标准，则表明该分子与周围氨基酸具有良好的结合状态。本发明发现，在28853蛋白质结构簇中，当考虑至少四个氨基酸残基与小分子接触时，小分子配体并没有位于合理的位置。此外，当研究单个簇成员时，在所有6310个簇中对98693链的53580个PDB ID蛋白上小分子结合口袋位于合理位置。这清楚地表明，在同一蛋白簇中的代表结构中识别所得小分子结合口袋，也同样适用于同一簇中其他结构链。除此之外，小分子与不同链结合的，共计314个蛋白簇、涵盖1984个PDB ID结构和3686个链，在本发明中不予以考虑。

序列比对，归一化和突变

利用ClustW程序，对同一聚类中不同蛋白链进行多个序列比对，这是排列同一簇中蛋白质链序列以识别蛋白质相似区域的重要方法。基于序列比对信息，将蛋白中位于同一位置、氨基酸不同的非保守性位点设定为单个序列中的突变位点。此外，利用DSSP程序定义了二级结构和几何特征的对蛋白质数据格式，用于校准5996个蛋白质簇和95007链的二级结构和氨基酸位点，以获得氨基酸序列和位置编码信息。氨基酸序列归一化：根据序列相似性保守程度和位点校准信息，对蛋白簇中所有链的氨基酸序列进行归一化和编码，以使所有蛋白质在相同位置具有相同的氨基酸编码。突变:结合氨基酸序列归一化编码信息的结构和几何特征，本发明可以获得每个蛋白质链相同的二级结构区域信息和氨基酸序列突变位点信息。

小分子结合口袋鉴定

在本发明中，基于几何构型，通过搜索与配体相邻氨基酸的方法进行小分子结合口袋的预测。为最大化囊括所有氨基酸，将每条链中距离生物活性配体分子重原子

以内的所有残基设定为分子结合袋。如图1所示，来自同一蛋白簇、不同的结构链中，形成结合袋的残基彼此之间存在显着差异。与蛋白4GTZ链A相比，由于小分子配体分子不同，4GTY蛋白A链中存在5个新的形成小分子配体的氨基酸(圆圈标记)。此外，如果同一链中不同配体上重原子的距离小于

则这些配体被视为一个整体的化合分子。而且，若同一蛋白质链上存在不同的结合口袋区域，并且口袋之间不同区域的重叠率高于>60％，这些不同结合位点被整合为同一个口袋。

图2和图3为相同聚类中不同蛋白结合口袋可视化。其中小分子重原子周围

范围内氨基酸重原子显示为球形；组成小分子结合口袋的氨基酸显示为棍棒型；(A/B)不同PDB ID编码晶体结构4GTZ和4GTY上A链小分子结合口袋图。

亚型结构簇

在前期“序列聚类”步骤中，从每个蛋白聚类中均挑选出一个模板链，从而使数据集更易于管理，并消除由于任何蛋白质的过度表达而可能进入分析的任何偏见。图2和图3显示整个绑定位点应分为子绑定位点类型，这意味着将Cluster1分为两个不同的子集群：Sub-cluster1A和Sub-cluster1B。如果代表性链和其他链之间的结合位点残基的重叠率>60％，则该簇中的结构仅包含一个整合的结合口袋。

结合精细化序列聚类和结合位点分析，本发明最终鉴定获取得到6693个亚型蛋白结构簇，包含135670个不同的蛋白结构，共计79244个蛋白质。

各亚型结构簇中模板结构挑选

对于各亚型结构簇，模板结构选时，优先级高的元素选择设定为单个模板结构：

拥有生物小分子配体>无突变>氨基酸序列长度最长>分辨率最高。

结构叠合获取配体位点整合体

通常使用具有相似拓扑结构的已知模板，在蛋白结构表面上探索“口袋”或“空腔”。基于模板的方法，通过全局拓扑搜索同源结构；而后，假定的结合位点可以在蛋白结构叠合后进行坐标转化。使用软件UCSF Chimera根据结合口袋主链重原子叠加同一簇中的所有蛋白质链，以获得蛋白质的叠加簇。如图4所示，每条链、相同位置处的结合位点氨基酸被收集为一个集合体，其被称为配体位点整合体。

图5为组成亚型聚类结构簇内配体结合口袋“整合体”周围氨基酸残基，利用UCSFChimera软件，基于骨架CA原子，将所有蛋白链叠合在模板结构中。不同蛋白链、相同编码位置上氨基酸以球棍型、白色显示。小分子以球棍型、不同元素显示为不同颜色。

结合位点氨基酸构象差异化分析

蛋白质的生物学功能，取决于氨基酸残基对于周围环境变化所表的柔性变化特征，包括小分子配体结合和化学修饰。更进一步讲，组成配体结合口袋上的重要氨基酸残基的柔性特征，需要对其进行准确和精确化的分析。

在各种氨基酸柔性特征研究中，各种详细信息均表明蛋白构象变化与RMSD波动之间存在高度正相关性。相同的正相关趋势也适用于B-factor值体系。温度因子B-factor受晶体优化提取方法、步骤和温度影响，精确的反映了蛋白真实的波动性，静态、动态和晶格紊乱程度。小分子配体的存在也会影响活性位点的温度因子B-factor值。非生物活性配体化学分子，也是蛋白晶体结构中扰动噪音的来源。突变位点总是小分子配体与结合位点之间的结合，同时对蛋白质骨架的稳定性产生影响。显而易见的是，包括ASA(溶剂可及表面)和空间位阻之内的几何特征，可以用作蛋白质柔氨基酸柔性的评判指标。旋转角定义为二面角的组合。在相同的结合位点集合中，从每条链的残基的不同旋转状态观察到的残基在整个过程中称为柔性，否则称为刚性。

综上所述，本发明共收集了6个表征因子，以综合性的鉴定结合位点残基的结构波动状况。最后，开发了一种数学统计公式来确定不同结构模型之间相同位置的结合位点残基的结构灵活性重要性。

RMSD and B-factor

RMSD：重原子位置均方根偏差(RMSD)值，是测量不同蛋白质体系上原子间平均距离的一种度量形式。在此，本发明计算了模板结构和和目标结构之间每个结合位点氨基酸上所有重原子、侧链和骨架主链RMSD值。包括ASP，GLU，PHE和TYR氨基酸在内的侧链末端对称原子被视为同一原子。

D-factor：在典型的蛋白质晶体分辨率参数中，高温度因子反映了重原子的低经验电子密度指数，反之亦然。温度因子值，被存储在各PDB格式文件内ATOM和HETATM坐标信息中，因此B因子可以反应氨基酸原子的真实静态或动态迁移率。但是，由于结晶过程中，应用了不同的结构优化程序，在不同已知结构中测量B因子的规模可能不同。因此，在本研究中，本发明使用标准化数据来代替原始数据，以比较不同蛋白质结构链温度B-factor因子之间差异。原子归一化的B因子值计算如下：

Average B-factor value

Sample standard deviation of B-factor value:

Each atom normalized B-factor value:

B_{normalized-single}＝(B-factor-B_average)/B_{normalized-all}

不同结构与模板结构之间K277氨基酸上的B因子和重原子均方根偏差RMSD差异较大为其侧链结构，骨架重原子差异较小。如图7和图8所示，氨基酸P305上所有重原子、K277的4骨架主链原子中仅观察到微小的结构波动。相较于模板结构，其它蛋白链上K277氨基酸上的侧链原子间的温度因子B-factor和Rmsd差异较为明显的，表明此处氨基酸较为柔性。

溶剂可及化表面积

溶剂可及化表面比率，是通过计算不同结构状态下、相同氨基酸位点的溶剂可及化表面积比值。如前文“接触比率和ASA比率”部分所述，利用DMS软件计算在每条链中形成配体结合位点的所有氨基酸的ASA值。在相同位置上的氨基酸若拥有更强的柔性特性，往往在溶剂可及表面积变化方面具有较大差异。如图9所示，蛋白C和蛋白A之间的差异较大，表明当氨基酸结构发生很大变化时，溶剂可及化表面积比率也会相应变化。

图5-12为配体结合口袋周围氨基酸结构差异化分析示意图，其中同一结构聚类上不同三个蛋白质，分别以图6的A、图7的B和图8的C命名。其中氨基酸K277和P305作为氨基酸位点，用来展示结构差异化表征因子RMSD(图7的B和图8的C)、B-factor、二面角(图7的B)、溶剂可及化表面积(图9的D)、突变信息(图10的E)和空间位阻(图11的图F)

空间位阻

对各配体分子与所有整体口袋中所有氨基酸之间的空间位阻进行计算和比较分析。设定各重原子之间接触参数，并利用UCSF Chimera软件计算各原子间空间位阻。若两重原子之间范德华半径之和减去它们之间的距离并减去可能存在的氢键对的余量，两重原子之间的结构被视为结构重叠或空间位阻:

Overlap_ij＝R_VDWi+R_VDWj–d_ij–allowance_ij

Overlap Cutoff＝1,H-bond Allowance＝0.6

为了检测配体分子与氨基酸之间空间位阻，采用阈值标准临界值为

浮动值为

(默认空间位阻参数标准分别为0.6和

)。较大的正截距值，表征空间位阻较为显著；而负截距值，表明两重原子之间并不存在空间位阻。不同结构上、相同位置的氨基酸之间若存在较大空间位阻，表明此处氨基酸构象变化大。如图11的F所示，随着残基K277的构象变化，蛋白质A链上的氨基酸K277与蛋白质C上的配体存在明显的空间位阻，表明K227氨基酸位点区域柔性变化显著。

突变体

突变引起的轻微结构扰动，会导致整个蛋白质发生变化。在同一蛋白簇中，通过多序列比对和基于序列比对获取的位点信息，鉴定各蛋白质链的突变位点信息。蛋白质序列中突变频率高，表明蛋白质在此处柔性较大。这些变化可能非常大，实际上超过了每对蛋白质之间差异的范围。

根据突变前后氨基酸的几何片段和理化性质的变化，将突变位点分为不同等级，并赋予不同的权重系数值。其中，突变前后氨基酸若具有完全相反的物理和化学特性，将对该位点赋予较高的权重系数，如图10的E和图13所示。为突变之前和之后具有完全相反的物理和化学特性的氨基酸分配更高的加权因子，如图10的E、图13-14所示。

二面角

氨基酸残基的二面角，是定义蛋白质复杂三维结构的重要构象性质之一。相同位置上、不同链上氨基酸之间的二面角的差异，反映了模板结构与对比结构之间的构象变化。二面角，主要用于表征小分子结合口袋周围氨基酸侧链基团的旋转状况。二面角差异越大，该部分的氨基酸的柔性就越强。在本发明的研究中，每个残基的二面角主要集中在侧链重原子上，并通过UCSF Chimera软件进行计算，如图12的G所示。

关键柔性氨基酸位点

在每一个蛋白结构聚类中，赋予可表征蛋白质氨基酸柔性的权重因子不同系数，公式如下：

ΔKey_score＝1.0*MS_{Steric clash}+1.0*MS_Mutations+0.9*MS_{Dihedral angles}+0.9*MS_RMSD+0.9*MS_ASA+0.8*MS_B-factor

其中，突变信息和蛋白-小分子间空间位阻，被认为氨基酸柔性的最主要位点，都被赋予1.0的重要系数。此外，反应不同氨基酸位相间差异的二面角、RMSD值、溶剂可及化表面积值，被赋予0.9的权重系数，B-factor值权重系数为0.8。最终，本发明根据不同权重系数表征因子计算公式，最后得到综合性的打分，并对组成配体结合位点的所有氨基酸进行打分和排名，排名靠前的氨基酸残基表明该位点结构位相差异较大。最终，本发明挑选排名靠前的5个氨基酸位点作为该蛋白质上配体结合口袋内的关键氨基酸柔性位点。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (12)

1.一种蛋白质小分子结合口袋上关键柔性氨基酸的识别方法，其特征在于：包括以下步骤：

步骤1、在蛋白质PDB结构中，通过95％的序列相似性将蛋白质的结构分类为不同的亚型簇；

步骤2、明确的生物学配体，将在同一结构亚型簇中所有蛋白结构重叠，并准确鉴定蛋白上小分子结合口袋；

步骤3、将已知3D结构的同一结构亚型簇中的相应结合位点，根据小分子结合位置划分为位点结合亚型；

步骤4、通过计算理化和结构参数差异表征因子来研究蛋白质柔性的性质；

步骤5、用形成蛋白质结合口袋的系统的最高评分来识别关键的柔性残基。

2.根据权利要求1所述的蛋白质小分子结合口袋上关键柔性氨基酸识别方法，其特征在于：

所述步骤1中，采用＞95％序列相似性阈值，利用Cd-hit软件对共计133976个PDB ID蛋白结构数据集，剔除冗余序列，从而得到总共41754个亚型结构簇和410079个非冗余链序列；

其中多蛋白链簇：具备多条蛋白链的集合簇，删除不符合条件的6591个数据簇集，最终得到35963个结构数据集。

3.根据权利要求1所述的蛋白质小分子结合口袋上关键柔性氨基酸的识别方法，其特征在于：

生物配体与蛋白质上邻近氨基酸残基的接触率，被视为配体与周围环境相互作用程度的几何量度；如果两个重原子的中心位于不超过范德华半径之和＞20％的距离处，则认为它们之间存在相互接触；

溶剂可及化区域(ASA)指小分子配体暴露于周围水溶液环境中的程度；DMS软件被用来计算41754个蛋白质簇中、不同结构状态下的生物活性化学分子的溶剂可及化区表面积值，公式如下：

ASA ratio＝1-ASA(complex)/ASA(apo)

如果小分子配体能够同时满足接触率≥60％和ASA率≤40％的阈值标准，则表明该分子与周围氨基酸具有良好的结合状态；在28853蛋白质结构簇中，当考虑至少四个氨基酸残基与小分子接触时，小分子配体并没有位于合理的位置；此外，当研究单个簇成员时，在所有6310个簇中对98693链的53580个PDB ID蛋白上小分子结合口袋位于合理位置。

4.根据权利要求1所述的蛋白质小分子结合口袋上关键柔性氨基酸的识别方法，其特征在于：

利用ClustW程序，对同一聚类中不同蛋白链进行多个序列比对；基于序列比对信息，将蛋白中位于同一位置、氨基酸不同的非保守性位点设定为单个序列中的突变位点；

利用DSSP程序定义了二级结构和几何特征的对蛋白质数据格式，用于校准5996个蛋白质簇和95007链的二级结构和氨基酸位点，以获得氨基酸序列和位置编码信息；

氨基酸序列归一化：根据序列相似性保守程度和位点校准信息，对蛋白簇中所有链的氨基酸序列进行归一化和编码，以使所有蛋白质在相同位置具有相同的氨基酸编码；

突变：结合氨基酸序列归一化编码信息的结构和几何特征，可以获得每个蛋白质链相同的二级结构区域信息和氨基酸序列突变位点信息。

5.根据权利要求1所述的蛋白质小分子结合口袋上关键柔性氨基酸的识别方法，其特征在于：

基于几何构型，通过搜索与配体相邻氨基酸的方法进行小分子结合口袋的预测；为最大化囊括所有氨基酸，将

来自同一蛋白簇、不同的结构链中，形成结合袋的残基彼此之间存在显着差异；与蛋白4GTZ链A相比，由于小分子配体分子不同，4GTY蛋白A链中存在5个新的形成小分子配体的氨基酸，以圆圈标记；此外，如果同一链中不同配体上重原子的距离小于

则这些配体被视为一个整体的化合分子；而且，若同一蛋白质链上存在不同的结合口袋区域，并且口袋之间不同区域的重叠率高于＞60％，这些不同结合位点被整合为同一个口袋。

6.根据权利要求1所述的蛋白质小分子结合口袋上关键柔性氨基酸的识别方法，其特征在于：

对于各亚型结构簇，模板结构选时，优先级高的元素选择设定为单个模板结构：

拥有生物小分子配体＞无突变＞氨基酸序列长度最长＞分辨率最高。

7.根据权利要求1所述的蛋白质小分子结合口袋上关键柔性氨基酸的识别方法，其特征在于：

使用具有相似拓扑结构的已知模板，在蛋白结构表面上探索口袋或空腔；基于模板的方法，通过全局拓扑搜索同源结构；而后，假定的结合位点可以在蛋白结构叠合后进行坐标转化；使用软件UCSF Chimera根据结合口袋主链重原子叠加同一簇中的所有蛋白质链，以获得蛋白质的叠加簇；每条链、相同位置处的结合位点氨基酸被收集为一个集合体，其被称为配体位点整合体。

8.根据权利要求1所述的蛋白质小分子结合口袋上关键柔性氨基酸识别方法，其特征在于：

RMSD：重原子位置均方根偏差(RMSD)值，是测量不同蛋白质体系上原子间平均距离的一种度量形式，在此计算了模板结构和和目标结构之间每个结合位点氨基酸上所有重原子、侧链和骨架主链RMSD值；包括ASP，GLU，PHE和TYR氨基酸在内的侧链末端对称原子被视为同一原子；

B-factor：在典型的蛋白质晶体分辨率参数中，高温度因子反映了重原子的低经验电子密度指数，反之亦然；温度因子值，被存储在各PDB格式文件内ATOM和HETATM坐标信息中，因此B因子可以反应氨基酸原子的真实静态或动态迁移率；但是，由于结晶过程中，应用了不同的结构优化程序，在不同已知结构中测量B因子的规模可能不同；因此，在本研究中，使用标准化数据来代替原始数据，以比较不同蛋白质结构链温度B-factor因子之间差异；原子归一化的B因子值计算如下：

Average B-factor value

Sample standard deviation of B-factor value：

Each atom normalized B-factor value：

B_{normalized-single}＝(B-factor-B_average)/B_{normalized-all}

不同结构与模板结构之间K277氨基酸上的B因子和重原子均方根偏差RMSD差异较大为其侧链结构，骨架重原子差异较小；氨基酸P305上所有重原子、K277的4骨架主链原子中仅观察到微小的结构波动；相较于模板结构，其它蛋白链上K277氨基酸上的侧链原子间的温度因子B-factor和Rmsd差异较为明显的，表明此处氨基酸较为柔性。

9.根据权利要求1所述的蛋白质小分子结合口袋上关键柔性氨基酸识别方法，其特征在于：

溶剂可及化表面比率，是通过计算不同结构状态下、相同氨基酸位点的溶剂可及化表面积比值，利用DMS软件计算在每条链中形成配体结合位点的所有氨基酸的ASA值。

10.根据权利要求1所述的蛋白质小分子结合口袋上关键柔性氨基酸识别方法，其特征在于：

计算和比较空间位阻：对各配体分子与所有整体口袋中所有氨基酸之间的空间位阻进行计算和比较分析；设定各重原子之间接触参数，并利用UCSF Chimera软件计算各原子间空间位阻；若两重原子之间范德华半径之和减去它们之间的距离并减去可能存在的氢键对的余量，两重原子之间的结构被视为结构重叠或空间位阻：

Overlap_ij＝R_VDWi+R_VDwj-d_ij-allowance_ij

Overlap Cutoff＝1，H-bond Allowance＝0.6

为了检测配体分子与氨基酸之间空间位阻，采用阈值标准临界值为

浮动值为

默认空间位阻参数标准分别为0.6和

较大的正截距值，表征空间位阻较为显著；而负截距值，表明两重原子之间并不存在空间位阻；不同结构上、相同位置的氨基酸之间若存在较大空间位阻，表明此处氨基酸构象变化大；随着残基K277的构象变化，蛋白质A链上的氨基酸K277与蛋白质C上的配体存在明显的空间位阻，表明K227氨基酸位点区域柔性变化显著；

鉴定突变体的突变位点：

突变引起的轻微结构扰动，会导致整个蛋白质发生变化；在同一蛋白簇中，通过多序列比对和基于序列比对获取的位点信息，鉴定各蛋白质链的突变位点信息；蛋白质序列中突变频率高，表明蛋白质在此处柔性较大；

根据突变前后氨基酸的几何片段和理化性质的变化，将突变位点分为不同等级，并赋予不同的权重系数值；其中，突变前后氨基酸若具有完全相反的物理和化学特性，将对该位点赋予较高的权重系数，为突变之前和之后具有完全相反的物理和化学特性的氨基酸分配更高的加权因子。

11.根据权利要求1所述的蛋白质小分子结合口袋上关键柔性氨基酸的识别方法，其特征在于：

计算二面角：每个残基的二面角主要集中在侧链重原子上，并通过UCSF Chimera软件进行计算；

确定关键柔性氨基酸位点

在每一个蛋白结构聚类中，赋予可表征蛋白质氨基酸柔性的权重因子不同系数，公式如下：

12.根据权利要求10所述的蛋白质小分子结合口袋上关键柔性氨基酸的识别方法，其特征在于：

突变信息和蛋白-小分子间空间位阻，被认为氨基酸柔性的最主要位点，都被赋予1.0的重要系数；

反应不同氨基酸位相间差异的二面角、RMSD值、溶剂可及化表面积值，被赋予0.9的权重系数，B-factor值权重系数为0.8；

根据不同权重系数表征因子计算公式，最后得到综合性的打分，并对组成配体结合位点的所有氨基酸进行打分和排名，排名靠前的氨基酸残基表明该位点结构位相差异较大。

Images