该文档翻译自Rosetta官网的教程性文档 Rosetta Tutorials, Demos, and Protocol Captures,原文链接:https://www.rosettacommons.org/demos/latest/Home 。“😇”内的内容为译者添加。
Contents
- 7. Scoring Tutorial
- 9. The Packer
- 9.1. The Problem of Optimizing Side-chains
- 9.2. ~ 9.3. Invoking the packer through the `fixbb` application & The sequence design problem
- 9.4. Working efficiently with the packer: TaskOperations
- 9.5. ~ 9.6. Protocols that use the packer & Calling and controlling the packer from RosettaScripts
- 9.7. How the packer algorithm works under the hood (for advanced users)
7. Scoring Tutorial
KEYWORDS: CORE_CONCEPTS ANALYSIS UTILITIES GENERAL
Tutorial by Shourya S. Roy Burman (ssrb@jhu.edu) Created 20 June 2016
Updated 29 May 2017 by Vikram K. Mulligan (vmullig@uw.edu) for the ref2015 energy function.
😇请初学者仔细阅读该教程;译者对该部分进行了结构上的重构,所以有些部分可能在顺序上与原文不匹配😇
Rosetta 使用基于物理和统计术语的能量函数来计算生物分子的能量。在本教程结束时,你应该理解:
- 在 Rosetta 中计算生物分子的能量意味着什么;
- 特定的相互作用类型对能量/得分的贡献有多大;
- 如何解释和比较 Rosetta 计算的能量分数;
- 如何制备用于评分的生物分子结构;
- 如何给一个生物分子打分;
- 如何将能量函数转换为其他预设函数;
- 如何出于某种目的定制一个能源函数中的术语;
- 如何得到每个残基对能量分数的贡献。
7.1. Score Function
在 Rosetta 中,生物分子的能量是通过评分来计算的。例如:Rosetta 有一个优化的能量函数或评分函数,称为ref2015
,用于计算由L-氨基酸组成的球状蛋白质中所有原子相互作用的能量;还有一些全原子评分函数用于其他生物分子,以及用于简化的质心表示的评分函数;此外,您可以创建一个自定义评分函数,以满足您的需求。
Rosetta 的评分函数是能量项的加权总和,其中一些代表物理力,如静电和范德华的相互作用,而另一些代表统计项,如在 Ramachandran 空间中找到扭转角的概率。以下是评分函数中使用的部分能量项列表(更多能量项解释见此处):
fa_atr 不同残基原子间的 Lennard-Jones 吸引力
fa_rep 不同残基原子间的 Lennard-Jones 斥力
fa_sol Lazaridis-Karplus 溶剂化能
fa_intra_sol_xover4 残基内部 Lazaridis-Karplus 溶剂化能
lk_ball_wtd 不对称溶剂化能
fa_intra_rep 同一残基中原子间的 Lennard-Jones 斥力
fa_elec 具有距离依赖性电介质的库仑静电势
pro_close 脯氨酸环闭合能量及其前残基 psi 角的能量
hbond_sr_bb 一级序列中相邻(相近)残基的骨架-骨架氢键
hbond_lr_bb 一级序列中较远残基的骨架-骨架氢键
hbond_bb_sc 侧链-骨架氢键能
hbond_sc 侧链-侧链氢键能
dslf_fa13 二硫键几何势能
rama_prepro Ramachandran 偏好(对于脯氨酸前位置和其他位置有单独的参数)
omega 骨架的ω二面角,为标准偏差为6度的平面性的调和约束。
p_aa_pp 给定的二面角为 phi 及 psi 的氨基酸的概率
fa_dun 由 Dunbrack 统计给出的侧链 rotamers 内能
yhh_planarity 一种使酪氨酸羟基保持在芳香环平面内的特殊扭势
ref 每种氨基酸的参考能量,平衡氨基酸项的内能,在设计环节发挥作用。
METHOD_WEIGHTS 本身不是能量项,而是每个氨基酸所使用 ref 能量项的参数。
ref2015
评分函数的具体权重值:
METHOD_WEIGHTS ref 0.773742 0.443793 -1.63002 -1.96094 0.61937 0.173326 0.388298 1.0806 -0.358574 0.761128 0.249477 -1.19118 -0.250485 -1.51717 -0.32436 0.165383 0.20134 0.979644 1.23413 0.162496
fa_atr 1
fa_rep 0.55
fa_sol 0.9375
fa_intra_rep 0.005
fa_elec 0.875
pro_close 1.25
hbond_sr_bb 1.17
hbond_lr_bb 1.17
hbond_bb_sc 1.17
hbond_sc 1.1
dslf_fa13 1.25
rama 0.25
omega 0.625
fa_dun 0.7
p_aa_pp 0.4
yhh_planarity 0.625
ref 1
7.2. Score Algorithm
Rosetta用来给一个结构打分的算法可以在这里找到。😇译者将这一部分翻译在了下面😇
EnergyMethod-related definitions
在 Rosetta 中,EnergyMethod 是打分的主要方法。一个 EnergyMethod 可以映射到多个打分项。例如,氢键的 EnergyMethod 对应的打分项有hbond_sr_bb
、hbond_lr_bb
、hbond_bb_sc
和hbond_sc
。
EnergyGraph
:存储相互作用的残基对的残基间能量;
EnergyEdge
:每条边都存储了单独一对残基的交互作用,由其构成EnergyGraph
;
EnergyMap
:将评分类型映射到它们的值;
ScoreFunction
:负责给一个姿态打分的主类。
One-body energies:残基内能量。它们与 EnergyGraph
分开存储。
Two-body energies:在残基对之间计算的残基间能量。它们存储在EnergyGraph
中。
Whole structure energies:整个结构的能量。
Context dependent:一个残基或一对残基的能量取决于周围残基的构象。例如,两个残基之间氢键的强度取决于每个残基相邻的数量。
Context independent:残基或残基对的能量不依赖于环境(例如 Lennard-Jones 能量)。
Short range:存在距离截断。超过这个距离,这个能量被认为是零(即,在两个残基之间能量不被评估,除非他们在这个距离截断内)。例如,氢键的能量范围很短。
Long range:没有距离截断。对所有适用的残基对进行能量评估。它们通常只在特定的残基对之间计算(例如 AtomPairConstraints
EnergyMethods Hierchy / EnergyMethods 层级 | ||
---|---|---|
One Body Energies | Context-independent one body energies | |
Context-dependent one body energies | ||
Two Body Energies | Short-range two body energies | Context independent short-range two body energies |
Context-dependent short-range two body energies | ||
Long-range two body energies | Context independent long-range two body energies | |
Context-dependent long-range two body energies | ||
Whole Structure Energies |
Overview of the scoring algorithm
- 确定哪些能量需要重新评分:对于相互移动的残基对,从
EnergyGraph
中删除这对残基的EnergyEdge
(那些能量需要重新评估);检测哪些残基对彼此接近(在 two body energies 的 short range 截止范围内);对于 a) 在距离截止范围内且 b) 它们之间的距离发生了变化的残基对,为这对残基添加一条EnergyEdge
到EnergyGraph
。 - 对于每一对残基,计算 short range 能量项:为所有残基对计算 context-dependent 的 two-body energies;如果残基相对于其他残基移动,则计算 context-independent 的 two-body energies;
- 对于每个残基:为所有残基计算 context-dependent 的 one-body energies;如果残基有内部自由度的变化(例如它的侧链被重新包装),则计算其 context-independent 的 one-body energies;
- 如果 long-range 的 two-body energies 是 context-dependent 且/或 残基发生相互移动,则对其重新打分;
- 为整体结构能量打分。
7.3. Comparing Rosetta Scores to Real-Life Energies
虽然Rosetta的大部分能量功能都是基于物理的,但它也有一些统计项以生成看起来像已知蛋白质结构的结构。虽然较低的评分结构更接近于原生结构,但评分并不能直接转换为像 kcal/mol 这样的物理能量单位。我们在此使用 Rosetta Energy Units (REU) 来代表得分的单位。此外,由于分数取决于所使用的分数函数,所以比较使用不同分数函数的结构的得分通常没有意义。
7.4. Demo
本教程中列出的所有演示命令都应该在下面目录下执行。这里所有的演示都使用linuxgccrelease
的安装版本。根据您的操作系统和编译器,您可能需要将其更改为其他相应的内容。
cd <path_to_Rosetta_directory>/demos/tutorials/scoring
7.4.1. Preparing PDBs for Scoring
😇该部分与教程 4.1 内容有较大重叠,可以略读;总的来说,Example 1. 回顾了 PDB 中存在 Rosetta 完全无法处理的分子的情况;Example 2. 回顾了PDB 中存在 Rosetta 可运行但将会重构的残基或分子的情况😇
Example 1.
要在 Rosetta 中对生物分子进行评分,我们使用score_jd2
可执行文件。与 Rosetta 中的大多数应用程序一样,该应用程序希望以某种方式格式化输入 PDB。从蛋白质数据库中直接下载的 PDB 可能适用于 Rosetta,也可能不适用于score_jd2
。这里有一个例子,我们试图给 PDB 3TDM 评分:
$ROSETTA3/bin/score_jd2.linuxgccrelease -in:file:s input_files/from_rcsb/3tdm.pdb
将会有报错:
ERROR: Unrecognized residue: PO4
这个 PDB 含有一个磷酸盐离子,如果没有额外的选项, Rosetta 无法处理。要对这个 PDB 进行评分,我们需要添加一个选项-ignore_unrecognized_res
以忽略 PDB 中的磷酸盐:
$ROSETTA3/bin/score_jd2.linuxgccrelease -in:file:s input_files/from_rcsb/3tdm.pdb -ignore_unrecognized_res
现在生成了一个score.sc
文件(此后后称为打分文件)。要进行下一步,请键入rm score.sc
。
Example 2.
如果一个输入 PDB 不符合 Rosetta 要求的确切规范,例如它缺少了 Rosetta 默认识别的重原子或非典型的残基(不像磷酸盐),Rosetta 会添加或改变原子以满足规范。您可以通过包含选项out:pdb
来要求它输出它实际得分的结构。在这个新的例子中,我们将对 PDB 1QYS 评分:
$ROSETTA3/bin/score_jd2.linuxgccrelease -in:file:s input_files/from_rcsb/1qys.pdb -out:pdb
在日志文件中会有:
...
core.io.pose_from_sfr.PoseFromSFRBuilder: Reading MSE as MET!
...
core.pack.pack_missing_sidechains: packing residue number 13 because of missing atom number 6 atom name CG
...
😇对log文件的解读略😇
需要注意的是:由于 Rosetta 会不确定地重建缺失的侧链,因此在 PDB 结构和评分文件方面,该示例每次运行都会产生不同的结果。
PDB 文件现在重建了缺失的原子并使用 MET 取代了 MSE,这是实际打分的结构。此外,打分文件将显示一个较大的正total_score
,表示一个不利的结构。这并不意味着结构是不稳定的,它只是意味着 Rosetta 相信在这个 PDB 中可能存在一些微小的空间冲突。分数文件的格式如下:
SEQUENCE:
SCORE: total_score dslf_fa13 fa_atr fa_dun fa_elec fa_intra_rep fa_rep fa_sol hbond_bb_sc hbond_lr_bb hbond_sc hbond_sr_bb linear_chainbreak omega overlap_chainbreak p_aa_pp pro_close rama ref yhh_planarity description
SCORE: 267.496 0.000 -422.275 290.201 -25.824 1.313 238.436 248.433 -1.045 -23.835 -2.245 -22.744 0.000 1.234 0.000 -4.258 0.000 2.749 -12.643 0.000 1qys_0001
为了避免这些问题,建议您总是使用与您最终打算使用的评分函数相同的relax
protocol来改进 PDB,这能减轻冲突,为 Rosetta 的成功做好准备。关于实现这一步的更多细节将在后面的教程中介绍。让我们把注意力转移到为经过优化的 PDB 评分上。要进行下一步,请键入rm score.sc
。
7.4.2. Basic Scoring & Analysis of the Score File
在本节中,我们将为 PDB 1QYS 评分(<path_to_Rosetta_directory>/demos/tutorials/scoring/input_files
提供了一个优化的版本)。首先,我们将使用默认的评分函数,即ref2015
。当我们开始使用更多的选项时,我们将开始使用 flags 文件(以后称为标签文件),而不是使用空格分隔选项列表。现在,我们在标签文件中输入的唯一选项是输入 PDB 和输出打分文件的名称:
-in:file:s input_files/1qys.pdb
-out:file:scorefile output_files/score.sc
运行下列指令:
$ROSETTA3/bin/score_jd2.linuxgccrelease @flag
将产生如下score.sc
:
SEQUENCE:
SCORE: total_score score dslf_fa13 fa_atr fa_dun fa_elec fa_intra_rep fa_rep fa_sol hbond_bb_sc hbond_lr_bb hbond_sc hbond_sr_bb linear_chainbreak omega overlap_chainbreak p_aa_pp pro_close rama ref time yhh_planarity description
SCORE: -163.023 -163.023 0.000 -423.638 109.662 -46.146 1.040 49.117 241.309 -3.934 -26.998 -11.234 -25.491 0.000 4.211 0.000 -13.603 0.000 -4.905 -12.643 1.000 0.230 1qys_0001
名为total_score
的第一列表示结构1QYS的加权总分。与前一节中的未细化结构相比,该结构的total_score
要低得多。对于这种规模的精细化结构,典型的评分为 -100 ~ -300 REU。分数越低,说明特定蛋白质的结构越稳定。在此有一个经验法则:当用ref2015
评分函数对一个精化结构进行评分时,每个残基得分 -1 ~ -3 REU 是常规的。
其他列代表计入总分的各个部分。例如,列fa_atr
表示不同残基原子间 Lennard-Jones 吸引力势能的加权分数。这种分解可以帮助确定哪些能量项比其他的贡献更大,即蛋白质中发生了什么样的相互作用。在这个具体的例子中,Rosetta 预测残基间的范德华力具有最大的稳定化贡献(-423.638),而溶剂化(以fa_sol
为代表)具有最大的去稳定化贡献 (241.309)。
当我们对未细化结构进行评分时,fa_rep
加权评分为238.436,而较大的fa_rep
加权分数(即比fa_atr
的稳定效果大得多)表明结构中存在冲突。但是在细化结构中fa_rep
为 49.117。约 200 REU 的减少表明弛豫步骤缓解了原始 PDB 中存在的小冲突。
7.4.3. More Scoring Options
1) Changing the Score Function
现在,我们将尝试使用一个评分函数来评分两个结构,该函数被优化为结合两个蛋白质,称为对接。我们还想将评分文件重命名为score_dock.sc
。你可以在<path_to_Rosetta_directory>/main/database/scoring/weights
中找到评分函数的详细列表。我们现在使用的 flags 文件内容为:
-in:file:l input_files/pdblist
-score:weights docking
-out:file:scorefile output_files/score_docking.sc
运行指令:
$ROSETTA3/bin/score_jd2.linuxgccrelease @flag_docking
得到的score_docking.sc
内容如下:
SEQUENCE:
SCORE: total_score score dslf_ca_dih dslf_cs_ang dslf_ss_dih dslf_ss_dst fa_atr fa_dun fa_elec fa_pair fa_rep fa_sol hbond_bb_sc hbond_lr_bb hbond_sc hbond_sr_bb linear_chainbreak overlap_chainbreak time description
SCORE: -89.596 -89.596 0.000 0.000 0.000 0.000 -143.190 5.640 -1.371 -2.577 3.929 62.290 -0.824 -5.653 -2.502 -5.338 0.000 0.000 0.000 1qys_0001
SCORE: -55.214 -55.214 0.000 0.000 0.000 0.000 -119.403 4.578 -1.412 -2.044 10.233 66.103 -3.048 -3.799 -3.360 -3.061 0.000 0.000 0.000 1ubq_0001
第一行(在标题行之后)给出 1QYS 的分数,第二行给出 1UBQ 的分数。在这个评分文件中,我们可以看到 1QYS 的total_score
与之前的运行相比发生了变化。但由于我们得到了相同的 PDB 文件,预测的实验稳定性应该是相同的。应当注意:比较不同评分函数产生的评分是没有意义的。
我们还可以看到,上例中的dslf_fa13
能量项没有了,取而代之的是描述二硫键几何结构的四个项。fa_atr
加权评分具有较小的绝对稳定贡献(尽管原始的、未加权的fa_atr
评分必须保持不变,因为它是相同的结构)。1UBQ 氨基酸数小于 1QYS,但总分较高,这是因为 Rosetta 计算每个残基的分数然后加和。这不应该被解释为 1UBQ 比 1QYS 更不稳定。应当注意:在不同蛋白质之间,一个结构的稳定性与总得分之间不存在很好的相关性。
2) Patch Files and Changing Term Weights
现在,假设我们想修改对接评分函数中一些项的权重。例如,我们可能希望降低一个更一般的评分项的权重,以支持更具体的一组评分项,或者为约束添加评分项。有三种方法可以实现:
- 创建一个自定义权重文件,并将路径传递给
-score:weights
- 添加补丁文件修改现有权重
- 从命令行设置特定项的权重
在下面的例子中,我们将演示后两个。我们使用补丁文件docking.wts_patch
,其内容为:
fa_atr *= 0.423
fa_rep *= 0.100
fa_sol *= 0.372
fa_intra_rep *= 0.000
fa_pair *= 0.000
fa_dun *= 0.064
hbond_lr_bb *= 0.245
hbond_sr_bb *= 0.245
hbond_bb_sc *= 0.245
hbond_sc *= 0.245
p_aa_pp *= 0.00
fa_elec = 0.026
dslf_ss_dst *= 1.0
dslf_cs_ang *= 1.0
dslf_ss_dih *= 1.0
dslf_ca_dih *= 1.0
pro_close *= 0.000
我们还将使用以下选项从命令行中将fa_atr
的权重重置为1.0(最初在对接时为0.338,通过docking.wts_patch
修改为 0.338 * 0.423 = 0.142974),在 flags 文件中:
-in:file:s input_files/pdblist
-score:weights docking
-score:patch docking
-score:set_weights fa_atr 1.0
-out:file:scorefile output_files/score_docking_patch.sc
运行:
$ROSETTA3/bin/score_jd2.linuxgccrelease @flag_docking_patch
得到的score_docking_patch.sc
内容为:
SEQUENCE:
SCORE: total_score score dslf_ca_dih dslf_cs_ang dslf_ss_dih dslf_ss_dst fa_atr fa_dun fa_elec fa_rep fa_sol hbond_bb_sc hbond_lr_bb hbond_sc hbond_sr_bb linear_chainbreak overlap_chainbreak time description
SCORE: -404.591 -404.591 0.000 0.000 0.000 0.000 -423.638 0.361 -1.371 0.393 23.172 -0.202 -1.385 -0.613 -1.308 0.000 0.000 1.000 1qys_0001
SCORE: -332.020 -332.020 0.000 0.000 0.000 0.000 -353.263 0.293 -1.412 1.023 24.590 -0.747 -0.931 -0.823 -0.750 0.000 0.000 0.000 1ubq_0001
注意,我们在补丁文件中将fa_pair
的权重设置为0。这消除了该项的任何贡献,因此打分文件缺少fa_pair
列。此外,根据补丁文件的定义,fa_rep
、fa_dun
等的加权分数也发生了变化。fa_atr
加权评分的增加是通过命令行对 flags 文件输入的权重增加的结果(顺便说一句,这与ref2015
评分函数的权重相同,因此fa_atr
的加权评分与最初的评分示例相同)。
7.4.4. Getting Individual Residue Scores
total_score本质上是Rosetta中各个剩余分数的总和。有时,对每个残差的评分进行详细分析可能更有用,这对于理解哪些残基有冲突特别有用。在本节中,我们将使用可执行文件per_residue_energies
和选项out:file:silent
来指定每个残基分解要写入的文件:
-in:file:s input_files/1qys.pdb
-out:file:silent output_files/per_res.sc
运行:
$ROSETTA3/bin/per_residue_energies.linuxgccrelease @flag_per_residue
得到的per_res.sc
内容如下:
SCORE: pose_id pdb_id fa_atr fa_rep fa_sol fa_intra_rep fa_elec pro_close hbond_sr_bb hbond_lr_bb hbond_bb_sc hbond_sc dslf_fa13 rama omega fa_dun p_aa_pp yhh_planarity ref score description
SCORE: input_files/1qys.pdb 3A -2.666 0.270 2.416 0.025 -0.269 0.000 0.000 0.000 0.000 -0.324 0.000 0.000 0.010 1.894 0.000 0.000 -1.630 -0.273 residue_1
SCORE: input_files/1qys.pdb 4A -5.618 0.237 2.802 0.026 -0.123 0.000 0.000 -0.564 0.000 0.000 0.000 -0.262 0.007 0.814 -0.348 0.000 1.081 -1.949 residue_2
...
第一行(在标题行之后)告诉我们链A中残基3的分数,在这里表示为3A(即列pdb_id
所示,PDB 残基 3A 是 Rosetta 编号下的残基 1,如residue_1
所示)。在文本编辑器中读取输出文件可能相当困难,因为列标题通常不直接位于列值的上方。你可能需要在 MS Excel, Libre Office Calc 或类似的电子表格程序中打开,以便更好地阅读它们。请确保在per_residue_energies
中使用与score_jd2
相同的评分函数,以获得相似的结果。
7.4.5. Getting Individual Residue Score Breakdowns
如果我们想进一步将这些分数分解为残基固有能量得分(one-body scores
)和相互作用得分,我们将使用residue_energy_breakdown
和以下选项:
-in:file:s input_files/1qys.pdb
-out:file:silent output_files/energy_breakdown.sc
运行:
$ROSETTA3/bin/residue_energy_breakdown.linuxgccrelease @flag_residue_energy_breakdown
得到的energy_breakdown.sc
内容如下:
SCORE: pose_id resi1 pdbid1 restype1 resi2 pdbid2 restype2 fa_atr fa_rep fa_sol fa_intra_rep fa_elec pro_close hbond_sr_bb hbond_lr_bb hbond_bb_sc hbond_sc dslf_fa13 rama omega fa_dun p_aa_pp yhh_planarity ref total description
SCORE: input_files/1qys.pdb 1 3A ASP -- -- onebody 0.000 0.000 0.000 0.025 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
...
SCORE: input_files/1qys.pdb 1 3A ASP 2 4A ILE -1.518 0.072 1.027 0.000 0.721 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.301 input_files/1qys.pdb_1_2
...
第一行(在标题行之后)告诉我们第一个残基(PDB 编号 3A)的内部能量分数(one-body scores
)的分解。后面还有一行表示残基 3A 和 4A 的相互作用能(如上所示)。同样,请确保在per_residue_energies
中使用与score_jd2
相同的评分函数,以获得相似的结果。
7.5. Tips
在今后的打分实践中,一定要牢记上文加粗的提醒部分。此外,有些蛋白质需要额外的文件来评分。例如,一个膜蛋白需要定义跨膜区域,而一个对称蛋白(使用对称框架)需要定义对称性。这样的 protocol 超出了本教程的范围。以下是这些方法的文档链接:
对 Rosetta scoring 的总结可以参考这篇文献:Alford et al. (2017) J. Chem. Theory Comput.
9. The Packer
KEYWORDS: CORE_CONCEPTS GENERAL STRUCTURE_PREDICTION DESIGN
Tutorial by Vikram K. Mulligan, Ph.D. Created 20 June 2016.
Rosetta优化侧链的主要算法被称为 Packer(封装器)。😇请初学者仔细阅读该教程;译者对该部分进行了结构上的重构,所以有些部分可能在顺序上与原文不匹配😇。在本教程结束时,你应该理解:
- 封装器解决的问题类型 ············ 😇 9.1.
- 如何使用
fixbb
应用程序调用打包程序。 - 如何控制封装器的性能。
- 如何确保封装器的最佳性能。
- 什么应用程序和算法调用封装器。
本教程还介绍了:
- 从RosettaScripts脚本语言调用封装器的常见方法。
- 使用RosettaScripts脚本语言控制封装器的常见方法。
- 封装算法如何工作的简要概述。
9.1. The Problem of Optimizing Side-chains
在 Rosetta 中,一个常见的任务是侧链的优化。我们可以这样思考这个问题:假设我们有一个结构 (我们将其称为姿态),其具有固定的主结构。在结构的每一个位置,我们都有一个关于侧链的离散可能性的列表,我们称之为 rotamers,其中一个 rotamer 是一个特定残基类型侧链的特定构象。我们想为每个位置选择一个 rotamer,这样的组合 rotamers 代表最低能量的解决方案。这就是封装器解决的问题。
但这个问题实际上是相当困难的:给定有 M 个残基的蛋白,其中每个位置的 N 种可能性,就有 NM 种可能。这很快就变成了一个天文数字的可能性——例如,3个 rotamer 在100个位置中的每个位置大约是 5x1047 种可能的组合,这无法穷尽枚举。为了解决这一问题,封装器使用了蒙特卡罗方法 (后面将详细讨论)。这意味着封装器是随机的;除非是最小的采样空间,否则它永远不会接近于穷尽式探索采样空间(的效果),并且返回的解虽然可能是一个好解,但不能保证是最好的解。需要注意的是:重复封装过程可能会产生一系列接近全局最优的类似解,这些都不一定是最好的解。
9.2. ~ 9.3. Invoking the packer through the fixbb
application & The sequence design problem
9.2.1. Basic Run
为了说明封装器的作用,让我们运行它。运行封装器的最简单方法是运行 Rosetta fixbb
应用程序:它存在的真谛(raison d’être 😇好高级的词汇😇)就是在输入姿态上调用封装器。在终端中,导航到demos/public/fixbb
目录,并运行如下命令(当然,你可能需要更改linuxgccrelease
以匹配你自己的版本;😇注意,这里使用了-resfile
标签,以对封装过程进行限制,关于这一点的解释可以在 9.3. 的结尾找到,但强烈建议你按顺序阅读,不必跳到 9.3. 提前了解😇):
cd <path_to_Rosetta_directory>/main/demos/public/fixbb
<path_to_Rosetta_directory>/main/source/bin/fixbb.default.linuxgccrelease -in:file:s 1l2y.pdb -in:file:fullatom -resfile resfile.txt -nstruct 5 >log.txt 2>err.txt &
这个应用程序包装了输入结构 (trp笼迷你蛋白,PDB 1L2Y) 的侧链。5次独立运行产生了5个输出结构。如果你将这些结构与输入结构进行比较 😇如上图,银白色为输入姿态,黄色及橙色为5个输出之二😇,你会发现 Rosetta 选择的侧链与输入结构相比略有不同;这是意料之中的,特别是考虑到 rotamer 的离散特性:真正“最好”的 rotamer 可能在两个被测试的 rotamer 之间,甚至可能永远不会被采样。
你可以选择完成演示的剩余部分 😇译者将这一部分作为 9.2.2. 翻译在了下面😇;这将告诉我们如何通过命令行和名为限制文件resfile的配置文件对封装器进行调整,以控制采样量、运行所花费的时间以及收敛到最优解决方案的可能性。
9.2.2. Application with Design
KEYWORDS: DESIGN GENERAL
Demo last modified by Vikram K. Mulligan, Ph.D. (vmullig@uw.edu) during the 2016 Documentation eXtreme Rosetta Workshop (XRW).
这是一个在固定骨架上运行的非常简单的设计 protocol 的演示。如果您以前从未运行过 Rosetta,那么这是一个可以运行的很好的第一个演示,因为它非常简单,选项很少。
😇Part 1😇
运行:
<path_to_Rosetta_directory>/main/source/bin/fixbb.default.linuxgccrelease -in:file:s 1l2y.pdb >log.txt &
运行后所得1l2y_0001.pdb
对比原姿态,可见 Rosetta 对序列及 rotamer 进行了替换,但骨架结构没有发生变化(银白色为原姿态,黄色为新生成的):
# 原序列
NLYIQWLKDGGPSSGRPPPS
EDKERHKKKGGONAGEPPQS
# 新生成的
😇Part 2😇
以序列配置文件的形式系统地列出序列变化;注意:必须将ROSETTA_TOOLS
环境变量设置为指向您的Rosetta/tools
目录,或者您可以手动输入Rosetta/tools
目录的位置;😇此外,似乎 Rosetta 的 Python 脚本都是用 2.x 版本写得,所以记得在一个 python2 的环境里运行😇:
ls 1l2y_0001.pdb > list.txt # 通常在一个列表中存放多个设计结构
python $ROSETTA_TOOLS/protein_tools/scripts/SequenceProfile.py -l list.txt -t 1l2y.pdb
然后要控制在每个序列位置允许哪些残基,您可以像这样添加一个限制文件(已包含在文件夹中);注意,我们将后缀变为_resout
,以不覆盖之前产生的文件。
rosetta/rosetta_source/bin/fixbb.default.linuxgccrelease -in:file:s 1l2y.pdb -resfile resfile.txt -out:suffix _resout > log_resout.txt &
9.3. The Sequence Design Problem
需要注意的是,如上所述的封装器问题没有对每个位置的候选侧链的性质作出假设;rotamer 可能性列表可以很容易地包含不同的侧链,就像它可以包含同一侧链的不同构象一样。因此,封装器是设计氨基酸序列的强大工具 (实际上,也是 Rosetta 的主要工具)。进入demos/public/fixbb_design
目录,并运行以下命令😇该命令与 9.2.2. 的第一步完全一致,结果也是相同的😇:
cd <path_to_Rosetta_directory>/main/demos/public/fixbb_design
<path_to_Rosetta_directory>/main/source/bin/fixbb.default.linuxgccrelease -in:file:s 1l2y.pdb >log.txt &
这将生成输出文件1l2y_0001.pdb
和score.sc
。如果你打开1l2y_0001.pdb
并将其与输入文件进行比较,您将看到序列发生了相当大的变化。然而,所选择的 rotamer 之间应该有合理的相互作用——也就是说,不应该有侧链占据相同的空间(碰撞)。注意,这里唯一的命令行选项是指定我们的输入文件;也就是说,在本例中,我们没有向fixbb
应用程序传递任何选项来控制封装器的行为。这就引出了非常重要的一点:封装器的默认行为是设计每个位置的侧链,这些侧链可以是全部20个标准氨基酸的所有 rotamer 之一。
fixbb_design
演示的剩余部分将教给你如何使用限制文件来控制包装器的行为,只允许设计特定的位置,并且只允许设计特定的氨基酸残基类型。😇该部分与 9.2.2. 完全相同,故不再翻译😇
😇在这里我们有必要就限制文件resfiles的内容进行一些讲解,在原教程中这是出现在 9.2.2. 的部分,但讲解的过于粗略。为了让读者理解为何 9.2.1. 中仅仅改变了 rotamer,而 9.2.2.第一步及9.3. 中连残基种类都改变了,我们首先来看一下resfile.txt
:
NATAA
start
其中:
标志 | 含义 |
---|---|
NATRO | 不处理该残基 |
NATAA | 保持在原位,但允许改变 rotamer |
ALLAA | 允许设计为其他任何氨基酸 |
PIKAA | 仅能设计为后面列出的氨基酸残基 |
例如 1 A PIKAA NT
代表允许将残基编号1改变为 N 或 T。
因此,使用了-resfile resfile.txt
标签的 9.2.1.教程中只改变了 rotamer;而在没有使用该标签的 9.3. 中不仅 rotamer 还有部分残基本身发生了改变😇
9.4. Working efficiently with the packer: TaskOperations
😇 这一部分可以不用细看,主旨就是一个合理利用限制文件 resfile;实际上后面的部分都不太需要细看了😇
当我们在没有设计的情况下运行fixbb
应用程序来重新包装侧链时,我们看到在速度、收敛概率和准确性之间存在制衡。包括更多的 rotamer 通常可以让 Rosetta 找到更好的侧链安排,但需要更长时间的运行且可能不太可能达到最低能量状态。封装器的设计大大增加了 rotamer 的数量:封装器必须考虑每一种侧链都有一个或多个 rotamer,而不是只包括一种侧链。在封装时,最好是进行短的测试,并注意在 tracer 输出中考虑的 rotamers 的数量。在 😇 9.2.1. basic run 😇 的例子中,氨基酸本身的类别是固定的,每个位置只考虑少量的构象,在输出的 log文件中有:
core.pack.pack_rotamers: built 256 rotamers at 20 positions.
由于 rotamer 数量少,位置少,封装工作几乎可以瞬间完成。后来,当我们打开额外的 rotamers (仍然保持固定的序列),输出信息改变为:
core.pack.pack_rotamers: built 6021 rotamers at 20 positions.
这仍然在Rosetta的可执行范围内,但此时,打包将需要一些时间(可能每次运行需要几秒)。没有额外 rotamers (这是我们在第二次演示中所做的)的设计也有一定的计算成本,但仍然是可控的:
core.pack.pack_rotamers: built 4737 rotamers at 20 positions
有多种不同的设置可以改变封装器的性能,从而使重构大量 rotamers 的工作更加高效。例如选项-linmem_ig 10
将使封装器在每个位置使用大量 rotamers 时更有效,但在每个位置使用较少数量时效率较低。然而,我们可以极大地简化问题,并通过resfile
使封装器在每个位置最多只有几个选择,就像我们在 9.2.2. 的第二部分中所做的那样:
core.pack.pack_rotamers: built 410 rotamers at 16 positions.
一般来说,为了保持搜索的快速和收敛性,尽可能地限制封隔器的选择是一个好主意。这就引出了另一个重要的概念:封隔器的每次下入都可以通过一个或多个TaskOperations
进行控制。
TaskOperations
可以通过几种方式之一传递给封装器,并且不同的TaskOperations
以不同的方式修改封装器的行为。我们已经看到了一些例子:用户可以通过命令行上的-resfile
选项指定限制文件来控制打包行为,它隐式调用了ReadResfile TaskOperation
。限制文件可以控制氨基酸种类、额外 rotamer 的存在以及由残基指数指定的基于 residue-by-residue 的封装器行为。有关限制文件及其完整特性的更多信息,请参见ReadResfile TaskOperation
文档和限制文件语法文档。其他的残基选项-ex1
、-ex2
等)也调用TaskOperations
来修改每个位置每个氨基酸类型的 rotamers 的数量。
TaskOperations
的一个重要属性是交换性:TaskOperations
可以以任何顺序应用,并且仍然产生相同的封装行为。为了实现这一点,某些由TaskOperation
控制的封装行为服从 AND 交换性:即如果TaskOperation
A、B和C允许该行为,那么当A、B和C同时应用时,该行为将被允许;如果A、B或C中的任何一个禁止该行为,那么当A、B和C同时应用时,该行为将被禁止。每个位置上允许的标准氨基酸类型服从 AND 交换性:例如,给定TaskOperations
A、B和C,只有当A、B和C均允许该位置为酪氨酸时,封装器才会将其设计为酪氨酸。换句话说就是:您只能禁止标准残基;且一旦禁止,它们就会保持禁止。
其他TaskOperation
功能服从 OR 交换性。例如,如果TaskOperations
A、B或者C打开了它们,就会出现允许额外的 rotamers。非标准残基类型服从 OR 交换性:当特定位置的特定标准残基类型的设计默认为开启,只能关闭时,特定位置的特定非标准残基类型的设计默认为关闭,只能开启。一旦允许一个非标准残基类型的,它就会一直允许。
最后,值得注意的是,TaskOperations
的互换性并不适用于用于控制最小化器的MoveMaps
(参见最小化器教程😇 将会在后续介绍😇 )。在MoveMap
中,后面的命令会覆盖前面的命令;在TaskOperations
的列表中,命令的效果是可交换的。
9.5. ~ 9.6. Protocols that use the packer & Calling and controlling the packer from RosettaScripts
封装器是一种基本的 Rosetta 算法,在许多更大的 protocol 中使用。许多 protocol 调用封隔器,包括abinitio
(在质心模式骨架构象搜索后使用封装器放置和优化全原子侧链)和relax
(在评分函数中增加排斥项的同时,交替进行一轮封隔和能量最小化),几乎每一种设计方案都使用封装器进行设计。
在RosettaScripts
中,可以直接使用PackRotamers
mover 调用封隔器。许多其他 mover,包括FastRelax
, FastDesign
和Disulfidize
也会调用封装器。通常,任何调用封装器的 Rosetta 组件都可以接收一个或多个TaskOperations
来控制封装器的行为。在RosettaScripts
中,TaskOperations
是在调用封装器的 mover 之前的脚本部分中单独声明的,然后通过名称传递给这些 mover。
9.7. How the packer algorithm works under the hood (for advanced users)
😇该部分为机翻,尚未检查😇
虽然使用Rosetta并不是严格要求详细了解封隔器的工作原理,但了解其中的一些原理可以帮助用户更有效地使用他或她的工具。要了解封隔器,首先要了解蒙特卡罗方法的概念。
广义地说,蒙特卡罗方法是一种试图用随机或随机移动来解决优化问题的方法。在实践中,罗塞塔的大多数蒙特卡罗方法都使用Metropolis- hastings算法,包括迭代步骤,在这些步骤中,一个人做出一个随机的移动,在某种程度上改变一个姿势,考虑一个姿势的能量变化作为这一移动的结果,然后根据称为Metropolis准则的规则接受或拒绝这一移动。Metropolis准则指出,如果移动导致能量减少,它总是被接受,而如果它导致能量增加,它被接受的概率等于e-ΔE/(kBT)。实际上,这意味着能量的增加越大,接受移动的可能性就越低,同时允许能量的小幅增加。
Metropolis标准包含一个“温度因素”,kBT,它决定了拒绝一个增加能量尺度的动作的概率如何随增加的幅度而增加。kBT的大值允许导致能量大幅增加的动作被频繁接受,而小值只允许导致能量小幅增加的动作被接受。虽然蒙特卡罗搜索可以使用固定的、任意选择的kBT值进行,但在模拟过程中改变这个值通常是有利的,在早期使用高值允许轨迹爬上障碍并逃离局部能量最小值,然后在稍后的模拟中降低它,以便轨迹“钻”到它所发现的最低能量井的底部。通常,许多kBT的上下循环可以让空间更完整的探索。涉及kBT项倾斜的蒙特卡罗方法广义上称为模拟退火方法。这种方法通常还包括有条件地后退到轨迹的早期部分,以逃避死胡同,以及其他增加收敛到最低能量状态概率的适应性行为。