【6.2.4】蛋白稳定性预测的工具比较（Provean/FoldX/ELASPIC/Rosetta）

protein_design 阅读量：次

预测突变对蛋白质的影响仍然是一个重要的问题。作为CAGI5 frataxin挑战的一部分，我们使用有限的八个突变数据集来评估Provean，FoldX和ELASPIC预测蛋白质的吉布斯自由能变化的准确性。我们发现不同的方法具有不同的优势和局限性，没有任何方法在所有指标上都严格优于其他方法。

  • ELASPIC达到了最高的准确性，同时还提供了可简化突变评估和分析的Web界面。
  • FoldX的准确性不如ELASPIC，但由于它不依赖于外部工具或数据集，因此更易于在本地运行。
  • Provean获得了合理的结果，同时在计算上比其他方法便宜，并且不需要蛋白质的结构。

除了提交给CAGI5社区实验的方法之外，为了高​​精度地告知其他方法，我们还评估了Rosetta的ddg_monomer协议，Rosetta的cartesian_ddg协议以及使用Amber软件包进行的热力学积分计算的预测。 ELASPIC仍然可以达到最高的精度，而Rosetta的catesian_ddg协议似乎在捕获数据的整体趋势方面表现最佳。

一、前言

DNA测序技术的进步导致可用基因组数据量的巨大增长。解释这些丰富的数据以提供有意义的和可行的见解仍然是一个挑战。基因组解释的一个重要方面涉及预测错义突变对蛋白质结构的影响。评估突变的结构影响以及蛋白质折叠的吉布斯自由能的相关变化（ΔΔG）可帮助预测突变的有害性（Glusman等人，2017，p.113），可提供一种机制解释特定突变如何产生特定表型（Nielsen et al。，2017），并有可能指导治疗策略的选择以及针对突变效应的靶向疗法的开发（Albanaz，Rodrigues，Pires，＆Ascher，2017）。尽管存在许多预测突变的ΔΔG的工具（Barlow等人，2018; Baugh等人，2016; Capriotti，Fariselli和Casadio，2005; Dehouck，Kwasigroch，Gillis和Rooman，2011; Kellogg，Leafer-Fay ，＆Baker，2011； Park等，2016； Pires，Ascher，＆Blundell，2014b； Schymkowitz等，2005），这些工具的准确性很难确定。大多数工具已经在通过实验测得的ΔΔG值的同一数据集上进行了训练和验证（Bava，Gromiha，Uedaira，Kitajima和Sarai，2004年），尽管它们通常报告该数据集具有良好的准确性，但结果却更多当涉及到以前未曾评估过的新突变时，其差异很大 (Buß, Rudat, & Ochsenreither, 2018; Geng, Xue, Roel‐Touris, & Bonvin, 2019; Khan & Vihinen, 2010; Kroncke et al., 2016; Potapov, Cohen, & Schreiber, 2009)。

基因组解释的关键评估（Critical Assessment of Genome Interpretation, CAGI）是一项社区实验，允许对与基因组变异的解释相关的任务进行客观评估和比较不同方法。在实验开始时，数据提供者会释放突变或整个基因组，并对其进行实验测量，并邀请来自世界各地的研究小组提交对这些实验测量的预测。在实验结束时，独立评估者将评估针对每个挑战的提交内容，并使用一组相关的指标来选择最准确的提交内容。 CAGI5 frataxin挑战的目标是针对人类frataxin（FXN）蛋白质的八个突变，预测吉布斯蛋白折叠的吉布斯自由能（ΔΔG）的变化。我们对此挑战提交了三份报告，其中包含Provean（Choi，Sims，Murphy，Miller和Chan，2012年），FoldX（Guerois，Nielsen和Serrano，2002年）和ELASPIC（Witvliet等人，2016年）的预测,其提交标识符分别为G6-3，G6-2和G6-1。由我们的实验室开发的计算框架ELASPIC做出的预测被独立评估者选为最准确的陈述。

在本文中，我们简要概述了预测突变ΔΔG的不同方法，并描述了我们对CAGI5 frataxin挑战提出的三项建议。我们还描述了通过三种替代方法进行的预测：Amber Thermodynamic Integration（Amber TI; Case等，2005），Rosetta’s ddg_monomer协议（Kellogg等，2011）和 Rosetta’s cartesian_ddg协议（Park等，2016）。尽管后者的预测没有在盲目挑战中进行评估，因此有更大的偏见潜力，但我们认为将这些预测包括在我们的分析中对于全面介绍最广泛使用的方法是必要的。

Provean是一种基于序列的方法，它使用氨基酸的保守模式来预测给定的突变是否可能有害。基于序列的方法的其他示例包括SIFT（Ng和Henikoff，2003），PolyPhen-2（Adzhubei，Jordan和Sunyaev，2001）和CADD（Kircher等人，2014）。虽然基于序列的方法不能直接预测突变的ΔΔG，但是基于序列的方法所做的预测通常与实验性ΔΔG的测量结果密切相关，因为有害的突变比良性突变更可能破坏蛋白质的结构稳定性（Berliner，Teyra ，Colak，Garcia Lopez和＆Kim，2014； Kroncke等，2016）。

FoldX和Rosetta是基于结构的工具的两个示例，这些工具使用统计势以及各种采样技术来预测突变的ΔΔG。统计潜力包括多种功能，包括使用分子力学力场计算出的能量，在高分辨率晶体结构中发现特定骨架构象和旋转异构体的概率，以及旨在改善实验值和预测值之间一致性的自定义例程的输出。这些功能以各种方式集成在一起以做出最终预测。现有的利用统计潜力的方法并没有尝试对蛋白质构象的大变化建模，假设蛋白质的骨架保持固定或仅允许突变位点周围的局部运动。

ELASPIC是由我们的实验室开发的元预测器(meta‐predictor)，它使用梯度增强决策树算法（Friedman，2002年）来整合Provean所做的预测，使用FoldX计算的经验能量项以及其他功能来预测ΔΔG突变。 ELASPIC属于使用序列和结构信息预测突变的ΔΔG的方法类别，该类别中其他方法的例子包括DUET（Pires，Ascher，＆Blundell，2014a），VIPUR（Baugh等人，2016）。 ）和STRUM（Quan，Lv，＆Zhang，2016）。在每种情况下，序列和结构特征都使用在实验测量的ΔΔG值的数据集上训练的机器学习算法进行集成（Bava等人，2004； Moal＆Fernández‐Recio，2012）。

Amber TI是一种协议，它使用Amber的pmemdGTI模块（Lee，Hu，Sherborne，Guo和York，2017）模拟从野生型到突变蛋白的转变并计算该转变的ΔΔG。它属于执行“炼金术”（alchemical）自由能计算的方法类别，在该方法中，分子动力学用于模拟从野生型到突变型蛋白质的跃迁，而跃迁的能量通过热力学积分（TI）多态Bennett计算接受率（mBAR，multistate Bennett Acceptance Ratio）或其他技术（Gapsys，Michielssens，Seeliger和de Groot，2015年; Lee等人，2017年）。使用炼金术自由能计算方法的主要优点是它们可以应用于更广泛的问题，例如预测突变对小分子结合的影响（Lee等人，2017）， d-氨基酸肽（Garton，Sayadi，＆Kim，2017）。炼金术自由能计算的主要缺点是，与其他方法相比，它们的计算量大得多。

虽然从只有八个突变的小型数据集得出结论时必须格外谨慎，但我们认为讨论属于这四个类别的每种方法的优势和局限性仍然是有益的。 需要更多的突变来进行盲目的评估，以证实或反驳我们的发现。

二、方法

2.1 FXN蛋白和突变

FXN是一种线粒体蛋白，参与铁硫（Fe-S）簇组装，血红素生物合成以及Fe2 +氧化为Fe3 +（Bridwell-Rabb，Winn，＆Barondeau，2011; O’Neill等，2005）。 FXN基因的破坏会导致弗里德赖希共济失调（Friedreich’s ataxia），这是一种神经和心脏变性疾病，具有常染色体隐性遗传方式，在欧洲人口中估计患病率为50,000（1）（Campuzano et al。，1996）。弗里德里希共济失调的大多数患者在FXN基因的第一个内含子中GAA重复扩增是纯合的，这会导致转录的FXN数量减少（Campuzano等，1996）。其余患者是GAA重复扩增和点突变的杂合子，已知有20多种不同的突变导致该疾病（Cossée等，1999； Galea等，2016； McCormack等，2000）。 GAA重复扩增和点突变杂合的患者可表现出多种表型，具体取决于存在的突变（Cossée等，1999）。对于错义突变，已经表明表型取决于突变对FXN的稳定性，FXN对铁硫装配复合物（SDU）的亲和力以及FXN的变构能力的影响。激活该复合物以促进铁硫团簇的生物合成（Bridwell-Rabb等，2011； Correia，Pastore，Adinolfi，Pastore和Gomes，2008）。

对于CAGI5 frataxin挑战，使用远紫外圆二色性和固有荧光光谱评估了人类FXN（NP_000135.2）中八个突变的热力学效应（Petrosino等人，2019）。 这些实验测量用于计算与突变相关的蛋白质折叠自由能的变化（ΔΔG）。 挑战的目标是使用计算工具来准确预测那些ΔΔG值。

2.2 ELASPIC Web服务器

Provean评分（G6-3），FoldXΔΔG值（G6-2）和ELASPICΔΔG值（G6-1）是通过将frataxin挑战的突变提交给ELASPIC Web服务器获得的（Witvliet等，2016）。有关挑战中所有突变的ELASPIC预测，请访问 http://elaspic.kimlab.org/result/a5393e/ 和表S1。 ELASPIC Web服务器使用FoldX版本3.0 beta 6.1，Provean版本1.1.5和ELASPIC版本0.1.42进行预测。对于更大的一组突变，也可以使用ELASPIC命令行实用程序获得预测（Strokach，Corbi-Verge，Teyra和Kim，2019年）。由于提供的实验值和FoldXΔΔG值是蛋白质展开而不是蛋白质折叠的吉布斯自由能的变化，因此这些值的符号相反。 Provean分数的符号也被颠倒了，因为Provean预测有害突变的值要比良性突变的值小，并且有害突变更有可能不稳定。

2.3 Rosetta

存在两种用于预测突变的ΔΔG的Rosetta协议：ddg_monomer协议（Kellogg等，2011）和Cartesian_ddg协议（Park等，2016）。 ddg_monomer协议使用具有soft_rep_design能量函数的灵活主干设计来生成50个野生型蛋白质模型和50个突变蛋白质模型。突变的ΔΔG计算为三个得分最高的野生型结构与三个得分最高的突变型结构之间的Rosetta能量之差（Kellogg等，2011）。 cartesian_ddg协议使用beta_nov16_cart能量函数在笛卡尔空间而不是扭转空间中优化了野生型和突变体结构。蛋白质的主链仅针对那些在6Å内或突变残基的三个氨基酸内的残基进行了优化。突变的ΔΔG计算为精制的突变体结构与精制的野生型结构之间的能量差，再乘以特定于能量函数的比例因子（Park et al。，2016）。使用Rosetta版本2017.26.59567生成ddg_monomer协议和cartesian_ddg协议的预测。表S2中列出了用于通过两种协议生成预测的系统命令。

2.4 Thermodynamic integration

热力学积分（TI）估计两个不同状态之间物理过程的自由能，其中哈密顿量H与参数λ链接，该参数λ用于将系统从状态A（λ= 0）转换为状态B（λ= 1； Seeliger＆de Groot，2010年）。所获得的自由能通过封闭的热力学循环的概念累加地组合起来，以获得计算出的自由能（Mitchell＆McCammon，1991）。折叠自由能差是通过未折叠状态模拟之间的自由能与为折叠蛋白质中的突变计算的自由能的组合而获得的。尽管蛋白质的未折叠状态很难建模，但已提出了各种方法来近似参考状态。在这里，通过建立一个模型来近似未折叠状态，在该模型中目标位置的侧翼残基被甘氨酸取代。这种方法被广泛用来模拟仿真过程中的展开状态（Seeliger＆de Groot，2010）。

预先平衡的结构倾向于生成更稳定的集合，因此可以更准确地估计自由能（Garton等，2018）。为此，使用PDB代码1EKG从结构中除去了所有水和离子原子。正确的质子化状态被识别并注释。在AMBER16（Case等，2005）和AMBER ff14SB力场（Maier等，2015）中使用TLEAP，通过添加一个12 nm3的显式水分子TIP3P来使结构溶剂化。接下来，添加Na +和Cl-抗衡离子以中和整个系统的净电荷，并应用周期性边界条件。然后，将它们最小化，平衡并加热至100 ps至300 K，并逐渐消除位置限制。用SHAKE约束氢键（Ryckaert，Ciccotti，＆Berendsen，1977），并使用2 fs的时间步长。粒子网格Ewald（Toukmaji，Sagui，Board和Darden，2000）算法用于处理远程相互作用。 5 ns后，完全消除了约束，并达到了完全平衡。然后，通过使用MMTSB工具集进行聚类来确定最具代表性的结构（Feig，Karanicolas和Brooks，2004年），用作后续TI的初始结构。

所有TI模拟都是在类似条件下进行的。除了一个步骤用于运动方程的积分外，仅在加热和预平衡步骤中停用了SHAKE算法。并且，使用9Å截止值进行的粒子网状Ewald方法（PME）进行了长距离静电相互作用。将λ= 0和λ= 1之间的转换分为11个窗口，其中λ值从0.0变为1.0，而Δλ= 0.1。用软核势处理整个突变残基，同时修改静电力和范德华力。首先在NVT集成中将所有起始结构最小化并在300 K时放松，这样可以节省物质的量（N），体积（V）和温度（T）。每个λ窗口的初始构象是通过每个λ值的1 ns预均衡顺序生成的。对于每个突变，每个λ窗口执行15纳秒的MD模拟。丢弃最初的1ns数据，并收集最后的14ns数据以500µfs的采样频率进行数据分析。每个模拟重复五次以计算集合平均值。在参考文献中可以找到有关推荐设置协议的更多信息（Garton等，2018; Lee等，2017; Seeliger＆de Groot，2010）。输入文件和脚本可以在这里找到： https://gitlab.com/kimlab/rapid 。

2.5 指标 Metrics

使用普通最小二乘法计算出图1中所示的最佳拟合线。 平均绝对误差是通过取期望值与实际值之间的绝对差的平均值来计算的。 使用SciPy的stats.pearsonr和stats.spearmanr函数计算Pearson和Spearman的相关系数。 使用scikit-learn的metrics.balanced_accuracy_score和metrics.roc_auc_score函数（Buitinck等，2013）来计算接收器操作员特征下的平衡精度和面积。 在计算平衡精度和接收器操作员特征下的面积时，将实验ΔΔG> 1 kcal / mol的突变归为不稳定，并将其分配为1，而实验ΔΔG≤1 kcal / mol的突变归为中性并分配 值为0。在计算平衡精度时，我们使用1kkcal / mol的阈值将突变分类为中性或不稳定。

三、结果

3.1 比较预测

图1显示了Provean（G6-3），FoldX（G6-2），ELASPIC（G6-1），Amber TI，Rosetta的ddg_monomer协议和Rosetta的cartesian_ddg协议的预测值与实验值之间的相关性，而图1中的残差最佳拟合线显示在图2中。Provean得分与实验ΔΔG值具有最强的相关性，皮尔逊相关系数为0.89，p值为0.003（图1）。高度相关的主要原因是Provean正确地捕获了p.Y123S和p.W173C突变的趋势，而其他方法则高估或低估了p.Y123S的影响，而低估了p.W173C的影响（图2）。突变p.Y123S的实验ΔG为4.48kcal / mol，与FXN的折叠ΔG接近，据报道为5.6kk / mol（Correia等人，2008）。然而，Provean预测p.Y123S的有害性不如p.W173C，这表明p.Y123S突变体仍比p.W173C突变体保留一些残留功能并且对生物体的危害较小。基于结构的方法难以为突变p.Y123S和p.W173C分配ΔΔG，因为它们无法模拟由那些突变引起的蛋白质构象的大变化。此外，大多数基于结构的方法都是使用实验测得的ΔΔG值的数据集进行训练或优化的，例如Protherm（Bava等人，2004），这些数据集仅包含ΔΔG高于5kkcal / mol的少数突变。 TI的琥珀色预测与实验ΔΔG值的相关性最差，皮尔逊相关系数为0.42，p值为0.3（图1d）。

除Rosetta的cartesian_ddg方案外，所有方法都高估了突变p.D104G（图2）的破坏稳定作用，该变化仅在ΔΔG为0.255kkcal / mol的情况下产生了轻微的破坏作用。突变p.D104G在α螺旋内部引入了一个甘氨酸，并且由于将与之相对应的phi-psi空间限制在较大的甘露糖占据的较小区域内，甘氨酸在α螺旋内部被大量消耗，甘氨酸可以被占用。 α螺旋结构的形成（Pace＆Scholtz，1998; Serrano，Neira，Sancho，＆Fersht，1992）。在p.D104G突变的情况下，由于其他结构特征，例如由从FXN的α螺旋突出的一组负残基引入的其他限制，此熵成本可能小于平均成本。现有的基于结构的工具将无法评估这些特征对α-螺旋形成的熵代价的影响，因为它们无法模拟蛋白质可以在折叠和展开状态下占据的完整构象集合。 Provean还预测，突变p.D104G的有害性要比纯粹基于其ΔΔG值所预期的有害。由于突变p.D104G从FXN表面的一小部分负残基中去除了天冬氨酸，因此可能会对FXN在将Fe2 +传递给血红素生物合成和从Fe2 +到Fe3 +氧化过程中的蛋白质中发挥的作用产生负面影响。这种获得的功能缺陷将与Provean的预测一致。

3.2评估方法性能

我们使用了几种不同的指标来量化由这六种方法得出的预测的准确性，并指出它们对于不同应用的适用性（图3）。方法的不同属性可能使它们比其他应用更适合某些应用。因此，我们在三种不同的子问题上评估了这六种方法的性能。首先，我们评估了这些方法能够准确预测突变引起的ΔΔG精确值的能力。在对热力学模型进行参数化或试图阐明突变导致特定表型的机制时，预测单个突变的ΔΔG可能很重要。其次，我们评估了该方法捕获一组突变的ΔΔG值的总体趋势的能力。捕获正确的数据趋势对于蛋白质设计等应用非常重要，在蛋白质设计中，经常使用计算方法来选择最可能稳定蛋白质或增加其与靶标亲和力的突变。最后，我们评估了该方法区分中性突变和不稳定突变的能力。如果我们要预测突变何时可能导致功能丧失或导致疾病，那么完成此任务的良好性能可能很重要。

参考资料

药企，独角兽，苏州。团队长期招人，感兴趣的都可以发邮件聊聊：tiehan@sina.cn
个人公众号，比较懒，很少更新，可以在上面提问题，如果回复不及时，可发邮件给我： tiehan@sina.cn