Rosetta Cartesian_ddG: 更快更准确的蛋白质稳定性预测方法

前言

计算单点突变自由能是分子稳定性改造的”老套路”。FoldX以及ddg_monomer application已经久经实战考验(具体用法参考前人总结:http://bioengx.hk1.91site.net/rosetta-ddg_monomer/)。今天我们就来说说Rosetta新一代的预测方法: Cartesian_ddG。

1 Cartesian_ddG方法概述

参考: Hahnbeom Park, Philip Bradley, Per Greisen Jr., Yuan Liu, Vikram Khipple Mulligan, David E Kim, David Baker, and Frank DiMaio (2016) “Simultaneous optimization of biomolecular energy function on features from small molecules and macromolecules”, JCTC.

Cartesian_ddG是一种新的采样方法，与ddg_monomer不同，对于骨架柔性的计算并不需要大量重复的强限制约束，而是采取了卡迪尔空间的优化来允许小幅度的骨架运动。

在计算∆∆G时，Cartesian_ddG方法使用两步Relax的方法:

-预先使用FastRelax进行优化，确定野生型以及点突变氨基酸最佳的侧链排布方式
-再使用卡迪尔空间版本的FastRelax对局部区域(点突变位置6埃范围内)进行彻底放松。
-评估野生型以及突变型蛋白的能量差，并使用能量函数特异性的校正因子对∆∆G值进一步矫正，以对应实验测定值单位kcal/mole。

效果：更快、更准。
Cartesian_ddG(使用opt-nov15函数)比ddg_monomer最佳Protocol(row16)要快上10倍有余。并且Cartesian_ddG预测∆∆G与实验∆∆G的皮尔森相关性从0.703升至0.743。

2 Cartesian_ddG的使用方法

Cartesian_ddG的使用方法十分简单。

-预优化野生型蛋白结构
-准备mutfile文件
-运行Cartesian_ddG
-分析∆∆G值

2.1 预优化野生型蛋白结构

创建relax_flag文件，文件内容如下:

-use_input_sc
-constrain_relax_to_start_coords
-ignore_unrecognized_res
-nstruct 20
-relax:coord_constrain_sidechains
-relax:cartesian-score:weights ref2015_cart # 务必设为ref2015_cart
-relax:min_type lbfgs_armijo_nonmonotone # 能量最小化的算法
-relax:script cart2.script
创建cart2.script文件，文件内容如下:(无需额外更改参数)

switch:cartesian
repeat 2
ramp_repack_min 0.02 0.01 1.0 50
ramp_repack_min 0.250 0.01 0.5 50
ramp_repack_min 0.550 0.01 0.0 100
ramp_repack_min 1 0.00001 0.0 200
accept_to_best
endrepeat
进行relax优化，打开终端并输入:

relax.mpi.macosclangrelease -s $pdb @relax_flag
打开relax的打分文件，查看最低能量的结构pdb编号，并以此作为出发的输入结构。

2.2 准备mutfile文件
单点突变的mutfile文件的格式如下:

注意！此处的89为pose序号，而非PDB中的残基号，设置错误一定会报错。
total 1 # mutfile中总突变数量
1 # 本轮突变氨基酸个数，单点突变设置为1.
T 89 A # 将89位的T突变为A
将文件保存为T89A.mutfile

2.3 运行Cartesian_ddG
Cartesian_ddG有两种运行模式: 蛋白稳定性预测模式(Protein stability mode)和 相互作用界面模式(Interface mode)。

1）蛋白稳定性预测模式(Protein stability mode)
创建cartddg_flag文件, 内容如下:

-ddg:iterations 5 # 默认为3，可以根据自身计算资源调整。
-ddg::cartesian
-ddg::dump_pdbs False # 是否输出突变后的蛋白PDB文件，推荐False,否则文件夹会很乱
-bbnbrs 1 # 骨架自由度，额外考虑线性邻居氨基酸的数目，1代表 邻近的3个氨基酸骨架自由度被考虑在计算中。
-fa_max_dis 9.0 # 控制范德华和溶剂化能量计算的范围，默认为6埃以内的氨基酸被考虑。
-score:weights ref2015_cart # 务必设为ref2015_cart
运行ddg计算，打开终端并输入，这么就会开始计算T89A的∆∆G。

cartesian_ddg.linuxgccrelease -s [inputpdb] @cartddg_flag -ddg:mut_file T89A.mutfile

2 ）相互作用界面模式(Interface mode)
该模式开发者并未做详细测试，斟酌使用，此处笔者未做测试。如有需要请参考: https://www.rosettacommons.org/docs/latest/cartesian-ddG

3 结果分析

运行完毕后，文件夹中会产生T89A.ddg这个文件。打开后，文件内容如下:

COMPLEX: Round1: WT: -333.782 fa_atr: -840.140 fa_rep: ………….
COMPLEX: Round2: WT: -333.782 fa_atr: -840.140 fa_rep: ………….
COMPLEX: Round3: WT: -333.782 fa_atr: -840.140 fa_rep: ………….
COMPLEX: Round4: WT: -333.782 fa_atr: -840.140 fa_rep: ………….
COMPLEX: Round5: WT: -333.782 fa_atr: -840.140 fa_rep: ………….
COMPLEX: Round1: MUT_89ALA: -330.352 fa_atr: -834.516 ………….
COMPLEX: Round2: MUT_89ALA: -330.352 fa_atr: -834.516 ………….
COMPLEX: Round3: MUT_89ALA: -330.352 fa_atr: -834.516 ………….
COMPLEX: Round4: MUT_89ALA: -330.352 fa_atr: -834.516 ………….
COMPLEX: Round5: MUT_89ALA: -330.352 fa_atr: -834.516 ………….
此处Round代表第几次迭代，根据-ddg:iterations设置定义。

WT:代表本轮能量计算，野生型的能量为-333.782。

MUT_89ALA:代表本轮能量计算，突变型的能量为-330.352。

我们需要将WT和MUT_89ALA平均后，用MUT能量平均值减去WT能量平均值就可以得到T89A的∆∆G值，为3.43。

稳定性变化评判标准:(务必根据自身需要调整)

∆∆G < -1.0 代表点突变为蛋白带来稳定效应； ∆∆G > 1.0 代表点突变为蛋白带来不稳定效应；

备注: 如果需要批量分析和点突变扫描，请使用脚本生成多个mutfile，然后按照以上流程进行计算。