导语：蛋白质在进化过程中，会发生氨基酸突变，破坏了残基之间的相互作用而导致蛋白质结构不稳定，若此时与突变残基具有相互作用的残基也随之发生突变，且构成新的相互作用，使蛋白质结构保持稳定，称这种变化为蛋白质共进化。本次介绍一篇由美国华盛顿大学霍华德休斯医学院研究团队于2014年5月发表在eLife上的一篇文章，该文章利用共进化信息实现了对蛋白质间残基-残基相互作用的鲁棒且准确的预测。
1
研究背景
蛋白质是连续的氨基酸序列经过盘曲折叠形成特定的三维结构。蛋白质必须有正确的形状才能正常运作，因为它们通常通过与其他蛋白质或分子结合而起作用。两个或多个蛋白质可以结合在一起，形成复合体来执行各种任务；研究这些复合物的结构是一项挑战，即使是在已知蛋白质亚基结构的情况下也是如此。作者注意到在蛋白质复合物中，来自不同蛋白质的氨基酸发生共进化时，这两个氨基酸往往会在蛋白质-蛋白质界面上发生接触。本文的作者开发了一种方法，利用共进化信息来预测在蛋白质复合物（two-protein complex）中哪一部分蛋白质发生了相互接触。
2
研究方法
(1)生成单个对齐
用HHblits和HHfilter对由EcoGene 3.0鉴定的4303个大肠杆菌蛋白基因进行多序列比对。为了减少冗余，作者为每个MSA（多序列比对）构造HMMs并且基于HHΔ（ 衡量HMM-HMM相似性：如果一对基因的HHΔ小于0.5，则划分为同一簇）聚簇基因，生成了2340个非冗余基因簇。
对于基准数据集， 使用与每个PDB关联的序列生成新的对齐。对于50S核糖体和NADH脱氢酶，作者使用来自PDB结构3uxr和4hea的嗜热菌HB8序列。
对于同源的NADH脱氢酶链L,M和N，作者在对齐生成协议中将e-value置为1E-60。除了来自大肠杆菌分析的复合物外，作者还使用来自PDB结构3ip4的序列将GatCAB酰胺基转移酶复合物纳入基准集中。对于PDB序列长度比平均长度范围长得多的情况，作者将覆盖范围过滤器修改为查询的50％。然后使用clustal omega v1.2重新比对序列。查询序列中不存在的残基不纳入后续分析中。
(2)生成配对对齐
作者从同一基因组中构建成对的蛋白序列[x1, x2, …, xp; xp+1, …, xp+q] ，位置1：p和p+1:p+q分别对应两个蛋白质。将这样一对蛋白质的多序列比对称为配对对齐。
对于基因组中有单拷贝基因的基因家族，如核糖体蛋白，很容易构建配对对齐，因为来自同一基因组的序列对可以直接连接。虽然在单基因组中一个基因存在的多个旁系同源基因的情况下，生成配对对齐的过程通常是复杂的，但在原核生物中，共同调控的基因通常在基因组中共同定位于操纵子中。作者限制使用具有小的，保守的基因间距离的基因对来创建配对序列，以此规避旁系同源基因。相似的方法被用于构建原核基因组中融合蛋白的数据库。将Δgene定义为一个基因对之间的注释基因数，我们仅考虑Δgene在60％的基因组中保守且小于20的基因对。考虑到在一个基因组中，大多数UniProt accession IDs是连续分配的，可以通过查看UniProt accession IDs的差异从而快速评估Δgene。然后对配对的比对进行过滤，以将冗余度降低至90％序列同一性，并去除缺口大于75％的位置。
(3)蛋白质复合体结构识别
为了在复合体结构中识别蛋白质对，对每一个大肠杆菌蛋白，用已生成的HHblits对齐构造HMM。接着使用hmmsearch在S2C数据库中扫描PDB序列。只考虑满足e-value小于1E-10的匹配。
(4)用配对对齐构建Gremlin模型
Gremlin为每个配对对齐构造一个全局统计模型，为配对对齐中的每一个氨基酸序列分配一个概率。
其中，Vi是编码位置特异性氨基酸倾向的向量，Wij表示编码位置i和j处氨基酸的耦合矩阵。通过最大化对齐的正则化伪似然度，从比对序列中获得这些参数，如下：
总和中的每一项是一个条件分布，该条件分布捕获了整个蛋白质序列中某个位置特定氨基酸的概率，R（v，w）是防止过度拟合的正则化项。
(5)用Gremlin评分排序残基对
为了将W ij矩阵简化为反映位置i和j之间耦合强度的单个值，作者首先计算S ij，即它们的向量2范数。使用这些值的行和列平均值来校正由于不同位置处的序列变异性而导致的Sij的差异：
其中，括号表示采用括号之外的指数的平均值，计算方法类似于Average Product Correction（APC）。与APC不同，作者仅通过计算对应于位置i和j的蛋白质位置的平均值来解决两个蛋白质家族中进化速率的差异：如果i和j都在第一个（第二个）蛋白质中，计算出第一个（第二个）蛋白质的位置上的平均值；如果i在第一个蛋白质中，而j在第二个蛋白质中，则仅在第一个蛋白质的位置上计算列平均值，而在第二个蛋白质的位置上仅计算行平均值。然后，计算归一化的耦合强度，ncsij,计算方法是用除以前3L/2的平均值。
作者观察到，当给定的复合物的最高耦合强度分数很高时，残基接触发生的频率越高越。为了解释这些依赖关系，建立了一个基于细菌50S核糖体复合物估计接触概率的模型：
通过与50S核糖体数据中观察到的频率的非线性拟合，确定了m、c和σ（分别为0.47、0.96和9.77）的值。
(6)将Gremlin评分转化成距离约束
作者将耦合强度转换为特定于残基对的距离约束，并将其包括在Rosetta结构预测程序中。作者使用以下形式的距离约束：
其中，d是受约束的原子间距，weight与 ncsij成比例。
(7)比较建模
使用RosettaCM基于与HHsearch生成的同源结构的比对（Remmert等，2011）建立了比较模型。对于预测会接触的区域中密度缺失的蛋白质，作者使用RosettaCM与受约束的共进化在对接之前构建缺失区域。
(8)从头建模
Rosetta从头建模的协议包括两个阶段：在初始阶段（“质心”），侧链由固定的质心原子表示，可以快速生成和评估各种蛋白质样拓扑结构；第二阶段（“全原子”）建立在明确的侧链中，并实现所有原子能的最小化。YIAM是一种膜蛋白， 由Rosetta膜能量函数建模。强排斥相互作用（公式1：weight: −100, cutoff: 35, slope: 2 and intercept: 100）被添加到细胞外区域的中心和预测的细胞内区域的中心之间，强吸引约束（weight:100, cutoff:35, slope:2 and intercept: 0）在预测的细胞内区域和细胞外区域内，有效地构建了膜状采样空间。使用MESSA的输出来预测跨膜区域。生成了100,000个模型，并且将最适合约束的20个模型收敛到单个群集。
(9)对接测试集
使用Jackhammer（HMMER v3.1b软件包的一部分）来识别基准集中18种复合物的子集，其中至少一种蛋白质或紧密同源物的apo形式具有可解析的结构。在结构为同源蛋白（e-value < 1E-20）并且存在大多数界面残基的情况下，我们使用比较模型生成了目标蛋白的结构模型。
(1)蛋白质-蛋白质对接
对于前3 / 2L预测中的每个约束间对，使用PatchDock v1.0，用聚类参数（rmsd 0.5; discardClustersSmaller 0）生成构象合集，然后使用所有约束对其进行评分。约束得分最高的5个模型使用约束在笛卡尔空间中进行了能量最小化。对于原生接触分数（Fnat）和界面均方根偏差（iRMSD）的计算，界面残基-残基接触是指任何重侧链原子之间的最小距离小于5Å。
3
实验结果
要识别两个蛋白质A和B之间的共进化残基对并不容易：只有当两个有机体中分别包含蛋白质A和蛋白质B的直系同源蛋白，并且生成A序列和B序列的对齐必须正确配对。为了简化直系同源鉴定，关注基因组中具有保守染色体位置的基因对，该基因对在基因组中被少于20个其他带注释的基因隔开。然后，作者为配对蛋白质家族中的序列建立Gremlin全局统计模型。接下来作者研究了具有较大耦合参数的残基对。
细菌50S核糖体亚单位中的残基-残基共进化
作者从研究细菌50S核糖体亚基中的残基-残基偶联参数开始，这是具有原子拆分结构的最大的进化保守细菌多蛋白复合物。对于复合物中的每个蛋白质，通过查询UniProt序列数据库构建多序列比对。对于每个这样的配对比对，建立一个Gremlin全局统计模型，计算归一化的偶联强度，并根据这些评分对蛋白质间残基对进行排名。偶联强度大于1表示两个残基之间的平均偶联率更高。
研究发现，在50S核糖体亚基中，只有一小部分残基协同进化，如偶联强度（图1A的 y轴）大于1.5时，并且大多数残基对距离在8Å 之内，所有残基对距离都在12Å之内。
图1A
图1B中显示了50S结构中共进化残基对的位置（为了清晰起见，各蛋白被拉开了）。黄线表示距离小于8Å，橙线表示距离小于12Å。对于50S核糖体，使用约1500个非冗余基因组的序列数据建立了Gremlin模型。
图1B
对于大的蛋白质-蛋白质复合物，复合物中蛋白质对之间的偶联强度的总和是否可用于区分直接相互作用的和非相互作用的蛋白质对？在50S亚基中，总偶联强度（图1C中的数字）大于1.5的蛋白质对，彼此相互作用（图1C中的方框）。但是，在50S亚基中有一些蛋白对接触，但没有发现共进化。显然，并非每种相互作用都可以通过偶联强度的总和来确定。
图1C
图1D表明，对于具有大量对齐序列的复合体，基于氨基酸序列共进化的残基-残基相互作用预测具有较高的置信度。
细菌复合体基准
对大肠杆菌的基因对生成配对对齐，对于1126个基因对，建立Gremlin全局统计模型，并确定了每个残基对的偶联强度。其中有64对基因对包含Gremlin得分>0.85的残基对，通过实验进一步确定了其中28对的三维结构，图2A中显示了其中一些复合物的Gremlin得分大于0.6的残基对的位置。
图2A
几乎所有Gremlin得分大于0.6的配对都在复合物结构中接触，除了NADH脱氢酶亚基（图2B）明显例外。有研究者认为该复合物在电子转移过程中经历了一系列构象变化。
图2B
结构未知的复合物的接触预测
在图3中，作者提供了目前结构未知的36种复合物的残基-残基接触预测。这些预测应该会有助于确定这些生物学上重要的复合物的结构。
图3
从接触预测到结构建模
预测的接触对于建模组建蛋白质复合物有帮助吗？作者在具有18种蛋白质复合物（复合物的部分结构已知）的对接测试集上进行评估。作者开发了一种对接协议，使用预测的接触作为距离约束，并采样了物理上合理的结构空间，以生成蛋白质-蛋白质复合物的模型。
图4A
iRMSD(界面均方根误差)最高的两种情况（图4A最后两行）在图4B和图4C中说明。高iRMSD是由于其中一种单体结构的构型在结合后发生了较大变化而引起的。尽管发生了这些变化，模型仍然可以准确识别结合界面。
图4B-C
4
总结
作者的研究结果表明，共同进化的残基对在蛋白质复合物中通常都会发生接触。当然，并不是所有在蛋白质界面发生接触的残基对都会发生共进化。之前的研究表明，只要有足够的比对序列，对单体蛋白的接触预测能达到很高的准确率，但是在这种情况下，必须先已知一个家族蛋白结构，并可以从中建立比较模型，限制了接触预测在结构预测中的效用。
作者使用提出的新方法对细菌中的28种蛋白质复合物中的蛋白质-蛋白质界面和细菌核糖体中蛋白质亚基之间的界面进行了预测，与真实结构比较后，得到了较高的准确率。下一步是考虑将该方法应用于真核生物中发现的蛋白质复合物。这是一个挑战，因为真核蛋白可供生成多序列比对的蛋白数量较少，因此更难检测到共进化现象的发生。
数据与工具
http://gremlin.bakerlab.org/complexes/
