RosettaVIP: 优化蛋白疏水内核折叠的稳定性设计方法
参考: https://www.rosettacommons.org/docs/latest/application_documentation/design/vip-app
教程资源: 位于$ROSETTA3/demos/public/vip 注意: demo中的结果格式已经过时，以本文实际测试结果为准。注意: demo中的flag比较诡异，以本文实际测试结果为准。
前言
RosettaVIP最早是为了解决人工设计蛋白不稳定的问题而开发的，RosettaVIP使用RosettaHoles的方法来发现蛋白质结构中折叠不合理的疏水空隙，并搜索合适的氨基酸来替代，使得蛋白核心折叠更加稳定。
演算法原理:
RosettaVIP设计是通过不断迭代以下3个步骤实现:
- 使用RosettaHoles搜索蛋白内部中的疏水空隙，并确定该空隙周围的氨基酸组成；
- 将这些氨基酸单点突变为疏水氨基酸(Fixed Backbone Design)；
- 将本轮具潜力的一些点突变进一步进行Relax优化的结构，计算单点突变的ddG值，并从这些点突变中选择最好的，在此突变的基础上进行下一轮迭代。
以上的迭代持续进行，直到找不到空隙，或没有更好点突变位点为止。
设计效果:
参考Automated selection of stabilizing mutations in designed and natural proteins一文
Benjamin 通过RosettaVIP演算法对数据集测试: 结果显示演算法相比于RosettaDesign演算法可以更加准确地识别蛋白疏水内核的稳定型残基。并以人工重新设计的λ受体蛋白DNA结合结构域、Protein L进行内核优化，试验验证结果提示RosettaVIP均能解决重设计蛋白不稳定性的问题。
此外，他们团队也尝试将RosettaVIP应用于天然蛋白的设计，以eMAP为例，通过设计获得了含5个点突变的eMAP5，其Tm值提升17.6℃。并且经过回复突变验证，这5个点突变只有组合在一起时，才能起到强大的组合效应。
1 准备输入文件
准备RosettaVIP的运行参数文件flags, 内容如下:
# 采样参数控制:
-s Start_Lambda_Model.pdb # 定义输入的结构
-ignore_unrecognized_res # 删除多余的无法识别的残基 以及 水分子
-cp:ncycles 4 # 控制迭代次数（不推荐使用）。
-cp:cutoff 6.0 # 定义可突变邻近氨基酸范围，文献值为7.0
-cp:max_failures 1 # 代表一次命令运行中，如果出现了x次错误，才会退出程序，设置为5次，可以确保RosettaHoles找到空隙。
-sasa_calculator_probe_radius 1.0 # 表面探针
# 控制侧链优化:
-cp:pack_sfxn ref2015 # 设置打分函数
-cp:exclude_file $filename # 允许定义不可突变的关键氨基酸
# relax优化控制
-cp:relax_sfxn ref2015 # 设置relax的打分函数
#-cp:skip_relax # 跳过relax阶段，只生成用于构建突变库的列表
-cp:relax_mover "relax" # 默认relax，可选"classic_relax"，但是更慢
#-cp:local_relax # 是否只进行局部relax的范围。针对大体系时启用，或对PDB结构进行截短处理，无关紧要的部分不考虑。demo中启用后，点突变寻找的能力下降，经常是失败的**。
# 报告输出
-cp:print_reports # 将生成 reports.txt报告文件。
-cp:print_intermediate_pdbs # 输出每一步的中间体蛋白PDB文件
exclude_file格式: 每行定义一个氨基酸(pdb number)。
- 如: 不允许A链128，136号计算突变:
128 A
136 A
写好之后，将文件的全名替换掉-cp:exclude_file后面的$filename即可。
2 运行RosettaVIP
demo以λ受体为例，测试RosettaVIP。
# 切换到工作目录
cd $ROSETTA/demos/public/vip/rosetta_inputs
# 直接运行
vip.mpi.macosclangrelease @flags
3 结果分析
运行结果输出至reports.txt文件:
可见，RosettaVIP对于我们输入的PDB结构总共进行了多轮迭代: 其中每轮都列出了候选的突变类型以及ddEgoe值(一种简化的能量函数打分值，可粗略地判断突变对稳定性的影响)
第一轮Relax后，程序分析出了10、13、35、52、57、60、66、71位点突变对结构稳定性有提升。而L13I的稳定性提升更明显(ddEgoe:-4.6042)，选取了该突变。
第二轮迭代筛选时是在L13I突变的基础上重新分析能够减少疏水空隙的点突变，新的结构基础上，突变的位置，以及能量贡献都发生了变化。
如此类推RosettaVIP模拟了"单点组合突变"的实验过程。
Iteration 1 : Found candidate mutations:
Position: 10 Native AA: ALA Mutant AA: LEU ddEgoe: -8.60701
Position: 13 Native AA: LEU Mutant AA: ILE ddEgoe: -1.75332
Position: 17 Native AA: PHE Mutant AA: ILE ddEgoe: -5.2701
Position: 31 Native AA: VAL Mutant AA: ILE ddEgoe: -1.65599
Position: 35 Native AA: ILE Mutant AA: LEU ddEgoe: -3.73465
Position: 42 Native AA: PHE Mutant AA: ILE ddEgoe: -7.19784
Position: 52 Native AA: PRO Mutant AA: LEU ddEgoe: -1.66479
Position: 57 Native AA: ALA Mutant AA: VAL ddEgoe: -5.48826
Position: 60 Native AA: PHE Mutant AA: ILE ddEgoe: -6.92573
Position: 61 Native AA: ALA Mutant AA: ILE ddEgoe: -5.89349
Position: 64 Native AA: PHE Mutant AA: MET ddEgoe: -9.8102
Position: 66 Native AA: VAL Mutant AA: ILE ddEgoe: -1.52301
Position: 68 Native AA: ILE Mutant AA: LEU ddEgoe: -1.13047
Position: 71 Native AA: PHE Mutant AA: TRP ddEgoe: -5.02144
Iteration 1 : The following mutations were accomodated after relaxation:
Position: 10 chain: 3 Native AA: ALA Mutant AA: LEU ddEgoe: -3.07001
Position: 13 chain: 3 Native AA: LEU Mutant AA: ILE ddEgoe: -4.6042
Position: 35 chain: 3 Native AA: ILE Mutant AA: LEU ddEgoe: -1.76918
Position: 52 chain: 3 Native AA: PRO Mutant AA: LEU ddEgoe: -1.53427
Position: 57 chain: 3 Native AA: ALA Mutant AA: VAL ddEgoe: -0.947129
Position: 60 chain: 3 Native AA: PHE Mutant AA: ILE ddEgoe: -2.34727
Position: 66 chain: 3 Native AA: VAL Mutant AA: ILE ddEgoe: -1.05786
Position: 71 chain: 3 Native AA: PHE Mutant AA: TRP ddEgoe: -3.56223
Accepted mutation from LEU to ILE at position 13 chain: 3
Iteration 2 : Found candidate mutations:
Position: 8 Native AA: ALA Mutant AA: ILE ddEgoe: -2.1497
Position: 10 Native AA: ALA Mutant AA: LEU ddEgoe: -8.0779
Position: 17 Native AA: PHE Mutant AA: ILE ddEgoe: -5.2137
Position: 31 Native AA: VAL Mutant AA: ILE ddEgoe: -1.65536
Position: 32 Native AA: ALA Mutant AA: MET ddEgoe: -4.9991
Position: 35 Native AA: ILE Mutant AA: LEU ddEgoe: -3.64756
Position: 42 Native AA: PHE Mutant AA: ILE ddEgoe: -7.00559
Position: 52 Native AA: PRO Mutant AA: LEU ddEgoe: -1.8697
Position: 57 Native AA: ALA Mutant AA: ILE ddEgoe: -5.61317
Position: 60 Native AA: PHE Mutant AA: ILE ddEgoe: -4.88678
Position: 61 Native AA: ALA Mutant AA: ILE ddEgoe: -5.65812
Position: 64 Native AA: PHE Mutant AA: MET ddEgoe: -6.68818
Position: 66 Native AA: VAL Mutant AA: ILE ddEgoe: -1.49476
Position: 68 Native AA: ILE Mutant AA: LEU ddEgoe: -1.04628
Position: 71 Native AA: PHE Mutant AA: TRP ddEgoe: -4.97966
.......................
.......................
对比最终效果的展示: 经过VIP design, 内部的疏水空隙明显减少。
特别注意事项:
- 不要额外调用
-ex1 -ex2选项，会大幅度下降计算速度。程序内部设计时已经考虑了额外的Rotamer； - RosettaVIP中relax是最耗时的部分，因此对大体系的计算能力不足，需要添加
-cp:local_relax选项。 - RosettaHoles1寻找疏水空隙是随机演算法，如果只允许一次，有可能找不到应有的空隙，因此需要多跑几次计算，可以尝试增加
-cp:max_failures的参数值。多次计算结果，会输出多条不一样的序列。 - 点突变的假阳率大约在25%左右。
备注: RosettaVIP在线伺服器:
https://rosie.graylab.jhu.edu/vip
推荐阅读：
