Notebook

MonteCarlo object & Movers¶

@Author: 吴炜坤

@email：weikun.wu@xtalpi.com/weikunwu@163.com

在上一个章节中，我们已经对Metropolis算法有了基本的了解。在本章节中，我们将介绍PyRosetta中的MonteCarlo对象以及介绍一些简单movers的使用。在最后我们将示例如何定义一个用于采样的方法流程。

一、MonteCarlo object¶

Metropolis算法有两个关键步骤，那就是move和accept，举例之前在mcmc.py中的代码函数, 判断接受过程中包括三个步骤, 记录move前的构象以及其能量, 采样后判断 $P_{accept}$ 是否接受新的构象。

在PyRosetta里，开发者们已经定义好了一个MonteCarlo的类，这个类主要有两个作用:

记录最低能量的Pose构象;
判断是否接受新的构象.

通过直接调用这个类和特定的moves组合使用，就可以实现Metropolis算法，以下我们将举一个实际的例子:

1 首先进行初始化和pose的生成:

In [1]:

# pyrosetta初始化
from pyrosetta import init, pose_from_sequence, create_score_function
from pyrosetta.rosetta.protocols.moves import MonteCarlo
init()
scorefxn = create_score_function('ref2015')

# 从序列生成三丙氨酸序列的线性Pose
pose = pose_from_sequence('AAA')

# 定义温度
kT = 1.0

# 定义MonteCarlo object:
mc = MonteCarlo(init_pose=pose, scorefxn=scorefxn, temperature=kT)

print(f'Old Pose Score:{scorefxn(pose)}')

PyRosetta-4 2021 [Rosetta PyRosetta4.conda.mac.cxx11thread.serialization.python37.Release 2021.26+release.b308454c455dd04f6824cc8b23e54bbb9be2cdd7 2021-07-02T13:01:54] retrieved from: http://www.pyrosetta.org
(C) Copyright Rosetta Commons Member Institutions. Created in JHU by Sergey Lyskov and PyRosetta Team.
core.init: {0} Checking for fconfig files in pwd and ./rosetta/flags
core.init: {0} Rosetta version: PyRosetta4.conda.mac.cxx11thread.serialization.python37.Release r288 2021.26+release.b308454c455 b308454c455dd04f6824cc8b23e54bbb9be2cdd7 http://www.pyrosetta.org 2021-07-02T13:01:54
core.init: {0} command: PyRosetta -ex1 -ex2aro -database /opt/miniconda3/lib/python3.7/site-packages/pyrosetta/database
basic.random.init_random_generator: {0} 'RNG device' seed mode, using '/dev/urandom', seed=2041883961 seed_offset=0 real_seed=2041883961 thread_index=0
basic.random.init_random_generator: {0} RandomGenerator:init: Normal mode, seed=2041883961 RG_type=mt19937
core.scoring.etable: {0} Starting energy table calculation
core.scoring.etable: {0} smooth_etable: changing atr/rep split to bottom of energy well
core.scoring.etable: {0} smooth_etable: spline smoothing lj etables (maxdis = 6)
core.scoring.etable: {0} smooth_etable: spline smoothing solvation etables (max_dis = 6)
core.scoring.etable: {0} Finished calculating energy tables.
basic.io.database: {0} Database file opened: scoring/score_functions/hbonds/ref2015_params/HBPoly1D.csv
basic.io.database: {0} Database file opened: scoring/score_functions/hbonds/ref2015_params/HBFadeIntervals.csv
basic.io.database: {0} Database file opened: scoring/score_functions/hbonds/ref2015_params/HBEval.csv
basic.io.database: {0} Database file opened: scoring/score_functions/hbonds/ref2015_params/DonStrength.csv
basic.io.database: {0} Database file opened: scoring/score_functions/hbonds/ref2015_params/AccStrength.csv
core.chemical.GlobalResidueTypeSet: {0} Finished initializing fa_standard residue type set.  Created 984 residue types
core.chemical.GlobalResidueTypeSet: {0} Total time to initialize 0.6941 seconds.
basic.io.database: {0} Database file opened: scoring/score_functions/rama/fd/all.ramaProb
basic.io.database: {0} Database file opened: scoring/score_functions/rama/fd/prepro.ramaProb
basic.io.database: {0} Database file opened: scoring/score_functions/omega/omega_ppdep.all.txt
basic.io.database: {0} Database file opened: scoring/score_functions/omega/omega_ppdep.gly.txt
basic.io.database: {0} Database file opened: scoring/score_functions/omega/omega_ppdep.pro.txt
basic.io.database: {0} Database file opened: scoring/score_functions/omega/omega_ppdep.valile.txt
basic.io.database: {0} Database file opened: scoring/score_functions/P_AA_pp/P_AA
basic.io.database: {0} Database file opened: scoring/score_functions/P_AA_pp/P_AA_n
core.scoring.P_AA: {0} shapovalov_lib::shap_p_aa_pp_smooth_level of 1( aka low_smooth ) got activated.
basic.io.database: {0} Database file opened: scoring/score_functions/P_AA_pp/shapovalov/10deg/kappa131/a20.prop
basic.io.database: {0} Database file opened: scoring/score_functions/elec_cp_reps.dat
core.scoring.elec.util: {0} Read 40 countpair representative atoms
core.pack.dunbrack.RotamerLibrary: {0} shapovalov_lib_fixes_enable option is true.
core.pack.dunbrack.RotamerLibrary: {0} shapovalov_lib::shap_dun10_smooth_level of 1( aka lowest_smooth ) got activated.
core.pack.dunbrack.RotamerLibrary: {0} Binary rotamer library selected: /opt/miniconda3/lib/python3.7/site-packages/pyrosetta/database/rotamer/shapovalov/StpDwn_0-0-0/Dunbrack10.lib.bin
core.pack.dunbrack.RotamerLibrary: {0} Using Dunbrack library binary file '/opt/miniconda3/lib/python3.7/site-packages/pyrosetta/database/rotamer/shapovalov/StpDwn_0-0-0/Dunbrack10.lib.bin'.
core.pack.dunbrack.RotamerLibrary: {0} Dunbrack 2010 library took 0.186065 seconds to load from binary
Old Pose Score:5.5084677688899575

2 尝试进行构象的move和接受判断:（重新运行代码之后，结构可能就已经变了）

In [2]:

# 对骨架二面角进行随机的扰动.
import random
new_phi = random.uniform(-180, 180)
new_psi = random.uniform(-180, 180)
pose.set_phi(seqpos=2, setting=new_phi)
pose.set_psi(seqpos=2, setting=new_psi)
print(f'New Pose Score:{scorefxn(pose)}')

New Pose Score:146.46348429332133

接下来尝试使用MonteCarlo对象进行接受判断:

In [3]:

# accept?
mc.boltzmann(pose)

Out[3]:

False

使用mc对象中的函数，可以给出更加详细的接受信息，如下，本次是第n次trials，接受的结果是？

In [4]:

# mc details;
mc.show_scores()
mc.show_counters()

protocols.moves.MonteCarlo: {0} MonteCarlo:: last_accepted_score,lowest_score: 5.50847 5.50847
protocols.moves.TrialCounter: {0}              unk trials=      1 NO ACCEPTS.

练习¶

请结合上述的函数，写一个随机采样骨架二面角的MonteCarlo采样程序(100次循环)。

二、Movers¶

A Mover is a type of object in Rosetta that interacts with a Pose. Frequently, a Mover changes the Pose.

从Mover的定义直观理解，所有能造成Pose中构象变化的操作，都可成为mover。在上述的例子中，我们通过单一地改变某一个二面角的操作，是蒙特卡洛的一次Move。接下来将会介绍几个与骨架优化相关的Mover，并且在下面的例子中读者将会学到如何将mover进行组装。

1 Small and Shear Moves¶

SmallMover和ShearMover是封装好的随机干扰phi/psi二面角的Mover，但有少许不同:

SmallMover: 随机选择一个氨基酸位点，扰动phi/psi二面角
ShearMover: 随机选择一个氨基酸位点，扰动phi角，随后将这个数值乘以-1，去干扰psi角。

SmallMover和ShearMover会同时进行骨架的合理性，确保干扰的残基位于Ramachandran plot允许的区域:

(图片来源: google图片)

这两个简单mover最常被用于小幅度扰动当前骨架结构的工具。

In [5]:

from pyrosetta.rosetta.protocols.simple_moves import ShearMover, SmallMover
from pyrosetta.rosetta.core.kinematics import MoveMap
movemap = MoveMap()
movemap.set_bb(True)
n_moves = 5  # 定义执行多少次随机扰动
kT = 2.0
SmallMover
small_mover = SmallMover(movemap_in=movemap, temperature_in=kT, nmoves_in=n_moves)
shear_mover = ShearMover(movemap_in=movemap, temperature_in=kT, nmoves_in=n_moves)

In [6]:

# small_mover/shear_mover 可以设定特定二级结构的扰动范围:
small_mover.angle_max(type="H", angle=25)
small_mover.angle_max(type="E", angle=25)
small_mover.angle_max(type="L", angle=25)

In [7]:

# 采样前pose的二面角分布:
print(pose.phi(1), pose.psi(1), pose.phi(2), pose.psi(2), pose.phi(3), pose.psi(3))

180.0 180.0 180.0 180.0 180.0 180.0

In [8]:

shear_mover.apply(pose)

core.scoring.ramachandran: {0} shapovalov_lib::shap_rama_smooth_level of 4( aka highest_smooth ) got activated.
basic.io.database: {0} Database file opened: scoring/score_functions/rama/shapovalov/kappa25/all.ramaProb
basic.io.database: {0} Database file opened: scoring/score_functions/rama/flat/avg_L_rama.dat
core.scoring.ramachandran: {0} Reading custom Ramachandran table from scoring/score_functions/rama/flat/avg_L_rama.dat.
basic.io.database: {0} Database file opened: scoring/score_functions/rama/flat/sym_all_rama.dat
core.scoring.ramachandran: {0} Reading custom Ramachandran table from scoring/score_functions/rama/flat/sym_all_rama.dat.
basic.io.database: {0} Database file opened: scoring/score_functions/rama/flat/sym_G_rama.dat
core.scoring.ramachandran: {0} Reading custom Ramachandran table from scoring/score_functions/rama/flat/sym_G_rama.dat.
basic.io.database: {0} Database file opened: scoring/score_functions/rama/flat/sym_P_rama.dat
core.scoring.ramachandran: {0} Reading custom Ramachandran table from scoring/score_functions/rama/flat/sym_P_rama.dat.
basic.io.database: {0} Database file opened: scoring/score_functions/rama/flat/avg_L_rama_str.dat
core.scoring.ramachandran: {0} Reading custom Ramachandran table from scoring/score_functions/rama/flat/avg_L_rama_str.dat.
basic.io.database: {0} Database file opened: scoring/score_functions/rama/flat/sym_all_rama_str.dat
core.scoring.ramachandran: {0} Reading custom Ramachandran table from scoring/score_functions/rama/flat/sym_all_rama_str.dat.
basic.io.database: {0} Database file opened: scoring/score_functions/rama/flat/sym_G_rama_str.dat
core.scoring.ramachandran: {0} Reading custom Ramachandran table from scoring/score_functions/rama/flat/sym_G_rama_str.dat.
basic.io.database: {0} Database file opened: scoring/score_functions/rama/flat/sym_P_rama_str.dat
core.scoring.ramachandran: {0} Reading custom Ramachandran table from scoring/score_functions/rama/flat/sym_P_rama_str.dat.

In [9]:

# 采样后pose的二面角分布:
print(pose.phi(1), pose.psi(1), pose.phi(2), pose.psi(2), pose.phi(3), pose.psi(3))
print(f'New Pose Score:{scorefxn(pose)}')

180.0 180.0 180.0 178.12064920489865 -178.12064920489865 180.0
New Pose Score:5.128540989020616

与之前的结果相比，可见使用SmallMover比我们随机在[-180, 180]的区间内随机选择一个数，构象的能量都要更低，构象也更加地合理，因为其检查了打分函数中的Rama项（骨架二面角势）

除了SmallMover和ShearMover以外，Rosetta的基础型Mover类型有非常的多，具体可参考官网链接: https://new.rosettacommons.org/docs/latest/scripting_documentation/RosettaScripts/Movers/Movers-RosettaScripts

2. MinMover¶

MinMover是Rosetta中被大量使用的一个Mover，基本的方法是使用梯度下降法优化。首先对Pose中每一个自由度做偏导计算梯度∇E, 然后以一定的步长去改变当前Pose的各个自由度分量，重新计算当前构象的能量。不断地迭代重复，直到能量收敛到忍受值范围之内。

举一个实际的例子说明: 当前体系中只有4个原子组成的二面角。通过4个原子的坐标可计算第1和第4个原子之间的范德华力得分。通过对两个原子距离的LJ势函数求导，可知两个原子之间的变化方向对能量的影响，选择梯度下降的方向对距离做一个很小的加量 $\delta d$ , 此时两个原子之间的距离为 $d+\delta d$ 。如此经过几轮迭代，当梯度∇E趋向于0且小于忍受值，迭代停止，构象达到能量较低的状态。

(图片来源: 晶泰科技团队)

MinMover是Rosetta中被大量使用的一个Mover通常与Metropolis Monte Carlo连用。MinMover有几个关键的设置，其使用的能量函数、定义Pose自由度的Movemap、梯度下降的方法以及能量的耐受值。

In [10]:

# 初始化定义MinMover
from pyrosetta.rosetta.protocols.minimization_packing import MinMover

# 定义movemap
my_movemap = MoveMap()
my_movemap.set_bb(True)

# 初始化minmover
min_mover = MinMover()
min_mover.movemap(my_movemap)
min_mover.min_type('lbfgs_armijo_nonmonotone')
min_mover.score_function(scorefxn)
min_mover.tolerance(0.01) # 能量变化的耐受值，当小于该值时停止优化.

关于min_type有比较多的选择: 可参考: https://new.rosettacommons.org/docs/latest/rosetta_basics/structural_concepts/minimization-overview

此处做出摘要和开发者的评论: "dfpmin" uses an exact line search, and Jim says it requires ~10 function evaluations per line search.

"dfpmin_armijo" uses an inexact line search, where the step along the search direction only needs to improve the energy by a certain amount, and also make the gradient a certain amount flatter (but not necessarily reach the minimum). This ends up being significantly more efficient, as once it gets going only 2-3 function evaluations are needed per line search, and approximately the same number of line searches are needed as for the exact dfpmin.

"dfpmin_armijo_nonmonotone" uses an even less exact line search along the descent direction, so that the step need only be better than one of the last few points visited. This allows temporary increases in energy, so that the search may escape shallow local minima and flat basins. Jim estimates this is several times more efficient than the exact dfpmin.

"dfpmin_strong_wolfe" uses the More-Thuente line minimization algorithm that enforces both the Armijo and Wolfe conditions. This gives a better parabolic approximation to a minimum and can run a little faster than armijo.

"lbfgs_armijo" uses the L-BFGS minimizer implementation with the Armijo inexact line search conditions.

"lbfgs_armijo_nonmonotone" uses the L-BFGS minimizer implementation with the even more inexact line search conditions.

"lbfgs_strong_wolfe" uses the L-BFGS minimizer implementation with the Wolfe conditions.

一般而言lbfgs_armijo_nonmonotone和dfpmin用的最多。

In [11]:

# 读入预先生成的高能量构象:
from pyrosetta.io import pose_from_pdb
bad_pose = pose_from_pdb('./data/HighEnergy.pdb')

core.import_pose.import_pose: {0} File './data/HighEnergy.pdb' automatically determined to be of type PDB

In [12]:

# 设置在Pymol中进行能量最小化的观察;
from pyrosetta.teaching import PyMOLMover
from pyrosetta.rosetta.protocols.moves import AddPyMOLObserver_to_conformation
pmm = PyMOLMover()
pmm.keep_history(True)
pmm.apply(bad_pose)
AddPyMOLObserver_to_conformation(bad_pose, True)

Out[12]:

<pyrosetta.rosetta.protocols.moves.PyMOLObserver at 0x7f8b88eea970>

In [13]:

# 进行能量最小化
min_mover.apply(bad_pose)

(图片来源: 晶泰科技团队)

观察可得能量最小化的轨迹。

3. Combination Mover¶

在Rosetta中有封装好的一些Combination Mover可以控制流程的运行逻辑，使用他们可以很方便的构建起一个MCMC流程。以下将着重介绍:

In [14]:

# 初始化
from pyrosetta.rosetta.protocols.moves import SequenceMover, RandomMover

SequenceMover和RandomMover都是可以通过add_mover来设置子Mover的。他们的控制逻辑有些不同:

SequenceMover: 按顺序将子Movers逐一执行;
RandomMover: 随机挑选一个子Movers执行;

In [15]:

# SequenceMover
seq_mover = SequenceMover()
seq_mover.add_mover(small_mover)
seq_mover.add_mover(shear_mover)
seq_mover.add_mover(min_mover)
seq_mover.apply(bad_pose)

In [16]:

# RandomMover
rand_mover = RandomMover()
rand_mover.add_mover(small_mover)
rand_mover.add_mover(shear_mover)
rand_mover.add_mover(min_mover)
rand_mover.apply(bad_pose)

TrialMover是将一个Mover(以下实例为small_mover)执行后与直接调用MonteCarlo对象，判断新构象是否被接受！

In [17]:

# demo for TrialMover #1
from pyrosetta.rosetta.protocols.moves import TrialMover

# 定义打分函数:
scorefxn = create_score_function('ref2015')
# 定义温度
kT = 1.0
# 定义MonteCarlo object:
mc = MonteCarlo(pose, scorefxn, kT)
trial_mover = TrialMover(small_mover, mc)
for i in range(5):
    trial_mover.apply(pose)
    mc.show_state()

protocols.moves.MonteCarlo: {0} MC: 1  3.59539  3.59539  3.59539  3.59539  0  0  0  accepted score beat low
protocols.moves.TrialCounter: {0}            Small trials=      1;  accepts= 1.0000;  energy_drop/trial=  -1.53315
protocols.moves.MonteCarlo: {0} MC: 1  3.66784  3.59539  3.66784  3.59539  0  0  0  accepted thermally
protocols.moves.TrialCounter: {0}            Small trials=      2;  accepts= 1.0000;  energy_drop/trial=  -0.73035
protocols.moves.MonteCarlo: {0} MC: 1  6.29155  3.59539  6.29155  3.59539  0  0  0  accepted thermally
protocols.moves.TrialCounter: {0}            Small trials=      3;  accepts= 1.0000;  energy_drop/trial=   0.38767
protocols.moves.MonteCarlo: {0} MC: 1  5.39026  3.59539  5.39026  3.59539  0  0  0  accepted score beat last
protocols.moves.TrialCounter: {0}            Small trials=      4;  accepts= 1.0000;  energy_drop/trial=   0.06543
protocols.moves.MonteCarlo: {0} MC: 1  5.42016  3.59539  5.42016  3.59539  0  0  0  accepted thermally
protocols.moves.TrialCounter: {0}            Small trials=      5;  accepts= 1.0000;  energy_drop/trial=   0.05832

它可以与SequenceMover或RandomMover组合来构建特定的mcmc采样逻辑:

In [18]:

# demo for TrialMover #2
# 随机选择min或small_mover+mcmc判断。
rand_mover = RandomMover()
rand_mover.add_mover(trial_mover)
rand_mover.add_mover(min_mover)
# 循环5次。
for i in range(5):
    rand_mover.apply(pose)
    mc.show_state()

protocols.moves.MonteCarlo: {0} MC: 1  5.42016  3.59539  5.42016  3.59539  0  0  0  accepted thermally
protocols.moves.TrialCounter: {0}            Small trials=      5;  accepts= 1.0000;  energy_drop/trial=   0.05832
protocols.moves.MonteCarlo: {0} MC: 1  5.42016  3.59539  5.42016  3.59539  0  0  0  accepted thermally
protocols.moves.TrialCounter: {0}            Small trials=      5;  accepts= 1.0000;  energy_drop/trial=   0.05832
protocols.moves.MonteCarlo: {0} MC: 1  5.42016  3.59539  5.42016  3.59539  0  0  0  accepted thermally
protocols.moves.TrialCounter: {0}            Small trials=      5;  accepts= 1.0000;  energy_drop/trial=   0.05832
protocols.moves.MonteCarlo: {0} MC: 1  5.42016  3.59539  5.42016  3.59539  0  0  0  accepted thermally
protocols.moves.TrialCounter: {0}            Small trials=      5;  accepts= 1.0000;  energy_drop/trial=   0.05832
protocols.moves.MonteCarlo: {0} MC: 1  5.42016  3.59539  5.42016  3.59539  0  0  0  accepted thermally
protocols.moves.TrialCounter: {0}            Small trials=      5;  accepts= 1.0000;  energy_drop/trial=   0.05832

同时也可以将已经定义好的RandomMover或SequenceMover作为基础mover输入TrialMover来构建更加复杂的逻辑:

In [19]:

# demo for TrialMover #3
# 先random，再mcmc判断。
top_mc = MonteCarlo(pose, scorefxn, kT)
top_trial_mover = TrialMover(rand_mover, top_mc)
for i in range(5):
    top_trial_mover.apply(pose)
    top_mc.show_state()

protocols.moves.MonteCarlo: {0} MC: 1  4.52295  4.52295  4.52295  4.52295  0  0  0  accepted score beat low
protocols.moves.TrialCounter: {0}        MoverBase trials=      1;  accepts= 1.0000;  energy_drop/trial=  -0.86365
protocols.moves.MonteCarlo: {0} MC: 1  4.5216  4.5216  4.5216  4.5216  0  0  0  accepted score beat low
protocols.moves.TrialCounter: {0}         MinMover trials=      1;  accepts= 1.0000;  energy_drop/trial=  -0.00135
protocols.moves.TrialCounter: {0}        MoverBase trials=      1;  accepts= 1.0000;  energy_drop/trial=  -0.86365
protocols.moves.MonteCarlo: {0} MC: 1  4.52058  4.52058  4.52058  4.52058  0  0  0  accepted score beat low
protocols.moves.TrialCounter: {0}         MinMover trials=      2;  accepts= 1.0000;  energy_drop/trial=  -0.00118
protocols.moves.TrialCounter: {0}        MoverBase trials=      1;  accepts= 1.0000;  energy_drop/trial=  -0.86365
protocols.moves.MonteCarlo: {0} MC: 1  2.38306  2.38306  2.38306  2.38306  0  0  0  accepted score beat low
protocols.moves.TrialCounter: {0}         MinMover trials=      2;  accepts= 1.0000;  energy_drop/trial=  -0.00118
protocols.moves.TrialCounter: {0}        MoverBase trials=      2;  accepts= 1.0000;  energy_drop/trial=  -1.50059
protocols.moves.MonteCarlo: {0} MC: 1  2.38253  2.38253  2.38253  2.38253  0  0  0  accepted score beat low
protocols.moves.TrialCounter: {0}         MinMover trials=      3;  accepts= 1.0000;  energy_drop/trial=  -0.00096
protocols.moves.TrialCounter: {0}        MoverBase trials=      2;  accepts= 1.0000;  energy_drop/trial=  -1.50059

上面的代码也可以改写使用RepeatMover来写

In [20]:

# demo for RepeatMover(顶层)
from pyrosetta.rosetta.protocols.moves import RepeatMover
top_mc = MonteCarlo(pose, scorefxn, kT)
top_trial_mover = TrialMover(rand_mover, top_mc)
rmover = RepeatMover(top_trial_mover, 5)
rmover.apply(pose)

特别值得注意的是，流程中不同层级的TrialMover可以独享或公用的一个MonteCarlo object。试想下，如果公用/独享一个MonteCarlo object会有什么差别？

综上所述，灵活使用这些Combination Mover可以构建出多种多样的搜索逻辑。

小练习¶

构建一个模拟退火的算法，算法流程如下:

内循环：

small moves 5次
Minimize
shear moves 5次
Minimize
mcmc判断

外循环：

重复5次，mc object的温度(kT)从2.5递减至0.5。

尝试使用上述的Mover，对一个20个氨基酸长度的螺旋序列进行骨架采样，重复1000-10000次，来预测蛋白质的结构。