蛋白质稳定性预测工具-Rosetta ddg_monomer

蛋白质由于其来自于生命体的本质，在工业环境下通常是不稳定的。暂不提高温高压或者含有有机溶剂的极端环境，就是在室温下抑或保存在4度的冰箱中，很多蛋白质都会很快的发生聚集、失活。

稳定蛋白质的方法有很多种，包括蛋白质工程，固定化，添加稳定剂等。蛋白质工程指的是对蛋白质进行基因改造，通过改变其结构来获得对包括pH，温度，有机溶剂等环境的耐受性。蛋白质工程通常可以采用三种策略，一种是定向进化，即在基因中随机引入突变，而后进行大规模的筛选，从而获得有益突变株。另一种策略是理性设计，即根据已知的蛋白质的结构与功能之间的关系，设计突变位点，而后通过定点突变的方式引入。最后一种策略，叫作计算机辅助设计。

顾名思义，计算机辅助设计主要依赖于计算机的计算。将蛋白质的三维结构作为输入，通过计算野生型与突变型的吉布斯自由能变的差值来判断突变型是否稳定。预测突变体是否稳定的软件有很多种[1]，包括CUPSAT、Dmutant、FoldX、I-Mutant、Eris [2]、Rosetta ddg_monomer [3]等。最近有研究者对这些方法的预测准确率进行了比较，发现Rosetta ddg_monomer，FoldX等的准确率要高于其他的软件[4]。

ddg_monomer已经成功的应用于包括柠檬烯环氧化物水解酶、卤代烷烃脱卤酶在内的多个蛋白质的稳定性工程化过程中。最近，我们也应用这个软件提高了大肠杆菌转酮醇酶的热稳定性，并对其预测准确率进行了评估，最后发现，其准确率可以达到65%[5]。面对着低效率的着定向进化，低成功率的理性设计，具有65%准确率的ddg_monomer 无疑是蛋白质工程学家的一个有力工具。由于本软件相关的中文资料并不多，BioEngX小编这里总结一下ddg_monomer的用法与操作步骤。

获取与安装

ddg_monomer程序是Rosetta软件内的一个应用。Rosetta软件最早是在华盛顿大学David Baker教授实验室开发的，目前软件内有多个应用可供用户使用，常用的应用程序有同源建模（comparative modelling）、短片段模拟与重建（Loop modelling and rebuilding）、蛋白质设计（protein design）、蛋白质与蛋白质对接（Protein-protein docking）、蛋白质配体对接（Protein-ligand docking）等。下载安装Rosetta软件后，这些程序就都可以应用了。

Rosetta软件对于学术界用户是免费的，只需要申请获得一个许可证，就可以从Rosetta的官网中下载软件了。详细的安装过程，大家可以从下面的链接中获悉：

https://www.rosettacommons.org/demos/latest/tutorials/install_build/install_build

需要注意的是，Rosetta软件目前只能安装在Linux系统或苹果Mac系统中。安装好Rosetta后可以在下面的目录中找到ddg_monomer应用程序rosetta/rosetta_source/src/apps/public/ddg/ddg_monomer.cc

输入文件准备

有两种方式可以用来运行ddg_monomer程序：一种是高分辨，一种是低分辨率。高分辨率的方式允许蛋白质的骨架（backbone）进行一定的运动，而低分辨率的方式只允许蛋白质侧链运动。根据文献的报道，高分辨率方式的准确率要比低分辨率方式高一些。本文主要介绍高分辨率方式运行ddg_monomer。

为了顺利以高分辨方式运行ddg_monomer 程序，我们需要准备以下文件：

• 对氨基酸进行重新排序的蛋白质PDB结构
• 距离限制文件
• Mutfile文件，该文件内用来指明需要突变的位点以及目标氨基酸

下面我们逐一介绍如何获得这三个文件以及每个文件的用途。

处理PDB文件

为了应用Rosetta软件，野生型蛋白的结构必须是已知的。从PDB数据中直接下载的PDB文件需要经过处理后，才能够应用于Rosetta 程序。处理的过程其实就是对PDB中的残基进行重新编号，将结构中第一个氨基酸编号为1。由于蛋白质内某些氨基酸灵活度比较高，比如N，C两端的氨基酸或某些无规则结构内氨基酸，在结构解析时这些残基经常会被漏掉，PDB文件内通常没有这些氨基酸的结构坐标。

Rosetta要求所有输入的结构文件，第一个残基都要被标识为残基1，而且整个结构内氨基酸排序不能有断点，所以大多数的PDB文件都需要进行处理，对氨基酸进行重新编号。Rosetta内提供了一个小的脚本程序，能够对PDB文件进行处理。可以通过下面的命令，运行这个程序

Python /usr/local/rosetta/tools/protein_tools/scripts/clean_pdb.py 1QGD AB > cleanpdb.log

上面的命令中，clean_pdb.py 是我们要调用的程序；1QGD是pdb文件的文件名；AB表示的是蛋白质的两条链，如果你的蛋白质中只有一条链，这里输入A就可以；cleanpdb.log 是本次程序运行的记录文件。

获得距离限制文件

高分辨率允许蛋白质骨架进行一定程度的运动。但为了防止骨架运动的幅度过大，严重偏离初始构象，需要对骨架的alpha C原子的位置进行一定的限制。距离限制也很简单，就是保证相邻的alpha C原子之间的距离在9 Å之间。Rosetta提供一个脚本，用来生成距离限制文件，可以采用下面的命令调用脚本：

tcsh /usr/local/rosetta/main/source/src/apps/public/ddg/convert_to_cst_file.sh mincst.log > input.cst

上面的命令中convert_to_cst_file.sh是我们的目标脚本，input.cst是本步骤生成的限制文件，mincst.log 是下一步能量最小化过程生成的log文件。

编写mutfile文件

Mutfile文件内需要标注清楚待突变氨基酸的位置以及目的氨基酸的种类。有关该文件的语法，感兴趣的同学可以从下面的网页中了解

https://www.rosettacommons.org/manuals/archive/rosetta3.4_user_guide/d1/d97/resfiles.html

小编在这里提供一个简单的例子。在一个txt的文件内输入下面的内容，而后将其扩展名更改为.mutfile，即获得了一个mutfile文件：

G 1 A

上面的内容表示，将PDB结构内的第一个氨基酸G突变成A。

能量最小化初始结构

野生型的三级结构需要进行能量最小化，而后才能用于ddg_monomer 的运算。Rosetta 提供了用于对PDB结构进行能量最小化的程序，其执行命令为：

/path/to/minimize_with_cst.linuxgccrelease -in:file:l lst -in:file:fullatom -ignore_unrecognized_res -fa_max_dis 9.0 -database /path/to/rosetta_database/ -ddg::harmonic_ca_tether 0.5 -score:weights standard -ddg::constraint_weight 1.0 -ddg::out_pdb_prefix min_cst_0.5 -ddg::sc_min_only false -score:patch rosetta_database/scoring/weights/score12.wts_patch > mincst.log

上面的命令可以对一系列PDB文件进行能量最小化，-in:file:l lst 用来指示需要进行能量最小化的文件，lst文件内标注待输入文件的文件名，可以是一个也可以是多个。这里需要强调的是，能量最小化输入的文件就是上一步处理过的PDB文件。Mincst.log文件是该命令的执行过程生成的log文件，这个文件可以通过脚本

convert_to_cst_file.sh

被转换为距离限制文件。

执行ddg_monomer程序

准备好输入文件后就可以正式执行ddg_monomer程序了。其执行命令如下：
/path/to/ddg_monomer.linuxgccrelease
-in:file:s WT.pdb
-ddg::mut_file H192P.mutfile
-ddg:weight_file soft_rep_design
-database /usr/local/rosetta/main/database/
-fa_max_dis 9.0
-ddg:minimization_scorefunction talaris2013
-ddg::iterations 50
-ddg::dump_pdbs true
-ignore_unrecognized_res
-ddg::local_opt_only false
-ddg::min_cst true
-constraints::cst_file input.cst
-ddg::suppress_checkpointing true
-in::file::fullatom -ddg::mean false
-ddg::min true
-ddg::sc_min_only false
-ddg::ramp_repulsive true
-unmute core.optimization.LineMinimizer
-ddg::output_silent

在上面的命令中，这三行指代输入文件：

-in:file:s WT.pdb 指的是能量最小化后的野生型的PDB结构，将WT换成你的文件名即可。

-ddg::mut_file H192P.mutfile 指的是mutfile文件，将H192P换成你的文件名即可。

constraints::cst_file input.cst 指的是距离限制文件的输入。

程序运行过程中，Rosetta软件会首先计算野生型的吉布斯自由能，而后计算突变型的吉布斯自由能，这个过程会持续50个循环。程序运行结束后，我们能够获得包括ddg_prediction.out 在内的一系列文件。ddg_prediction.out这个文件是最重要的，因为里面有我们需要的ΔΔG数据，也就是野生型吉布斯自由能与突变性吉布斯自由能的差值。如果ΔΔG数值为负，即突变性的自由能低于野生型的能量，那么说明Rosetta软件预测该突变型是有效的，能够提高野生型的稳定性。

上面我们介绍了如何运行ddg_monomer来预测单个点突变的稳定性。尽管小编认为整个过程已经介绍的很详细了，但在真正的执行过程中，通常还会遇到这样那样的问题。遇到问题怎么办呢？首先可以查看Rosetta提供的用户指南，查看这个链接https://www.rosettacommons.org/manuals/archive/rosetta3.4_user_guide/d3/d28/ddg_monomer_application.html，或参考文献[3]。

其次就是加入到我们BioEngX的大家庭，与其他同行交流讨论。感兴趣的小伙伴可以添加管理员微信号bioengxadmin，加入到BioEngX-分子动力学模拟 讨论群中。

Reference

[1] Khan S, Vihinen M. Performance of protein stability predictors[J]. Human mutation, 2010, 31(6): 675-684.

[2] Yin S, Ding F, Dokholyan N V. Modeling backbone flexibility improves protein stability estimation[J]. Structure, 2007, 15(12): 1567-1576.

[3] Kellogg E H, Leaver‐Fay A, Baker D. Role of conformational sampling in computing mutation‐induced changes in protein structure and stability[J]. Proteins: Structure, Function, and Bioinformatics, 2011, 79(3): 830-838.

[4] Thiltgen G, Goldstein R A. Assessing predictors of changes in protein stability upon mutation using self-consistency[J]. PloS one, 2012, 7(10): e46084.

[5] Yu H, Yan Y, Zhang C, et al. Two strategies to engineer flexible loops for improved enzyme thermostability[J]. Scientific Reports, 2017, 7.