Notebook

Fragment & Folding¶

@Author: 吴炜坤

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

经过前面的学习，读者应该已经掌握了使用mcmc算法进行蛋白质构象搜索的基本方法。此处我们尝试对一段螺旋序列进行结构的预测，还是使用我们之前的ShearMover/SmallMover/MinMover以及各种CombineMover来构建一个蛋白质Folding的程序:

一、崎岖的能量面¶

首先从简单的开始，我们使用PDB:ID为2i9m的一段alpha螺旋结构:

In [1]:

# load
from pyrosetta import init, pose_from_pdb

init()
pose = pose_from_pdb('./data/2i9m.pdb')

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=1287931986 seed_offset=0 real_seed=1287931986 thread_index=0
basic.random.init_random_generator: {0} RandomGenerator:init: Normal mode, seed=1287931986 RG_type=mt19937
core.chemical.GlobalResidueTypeSet: {0} Finished initializing fa_standard residue type set.  Created 984 residue types
core.chemical.GlobalResidueTypeSet: {0} Total time to initialize 0.637641 seconds.
core.import_pose.import_pose: {0} File './data/2i9m.pdb' automatically determined to be of type PDB
core.io.pose_from_sfr.PoseFromSFRBuilder: {0} [ WARNING ] discarding 1 atoms at position 1 in file ./data/2i9m.pdb. Best match rsd_type:  SER:NtermProteinFull
core.conformation.Conformation: {0} [ WARNING ] missing heavyatom:  OXT on residue GLY:CtermProteinFull 17

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

In [2]:

# 获取它的序列;
sequence = pose.sequence()
print(sequence)

# 读入预先准备好的线性多肽;
linear_pose = pose_from_pdb('./data/linear_pep.pdb')

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

In [3]:

# 使用TrailMover和其他Mover一起构建一个folding的程序;
# pyrosetta初始化
from pyrosetta import create_score_function
from pyrosetta.rosetta.protocols.moves import MonteCarlo
# 创建全原子打分函数:
scorefxn = create_score_function('ref2015')
# 定义温度
kT = 1.0
# 定义MonteCarlo object:
mc = MonteCarlo(linear_pose, scorefxn, kT)

# 定义movers
from pyrosetta.rosetta.protocols.simple_moves import ShearMover, SmallMover
from pyrosetta.rosetta.protocols.minimization_packing import MinMover
from pyrosetta.rosetta.core.kinematics import MoveMap
movemap = MoveMap()
movemap.set_bb(True)
n_moves = 1  # 定义执行多少次随机扰动
small_mover = SmallMover(movemap, kT, n_moves)
shear_mover = ShearMover(movemap, kT, n_moves)
small_mover.angle_max(25)
shear_mover.angle_max(25)

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

# 初始化combine mover
from pyrosetta.rosetta.protocols.moves import RandomMover
# RandomMover
rand_mover = RandomMover()
rand_mover.add_mover(small_mover)
rand_mover.add_mover(shear_mover)
# rand_mover.add_mover(min_mover)

from pyrosetta.rosetta.protocols.moves import TrialMover
trial_mover = TrialMover(rand_mover, mc)

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
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.155837 seconds to load from binary

好！上面已经准备好了mcmc的程序了，读者可以尝试运行10、100、1000次。查看最后的构象状态。

In [4]:

# 记录一下构象变化轨迹:
from pyrosetta.teaching import PyMOLMover
from pyrosetta.rosetta.protocols.moves import AddPyMOLObserver_to_conformation
pmm = PyMOLMover()
pmm.keep_history(True)
pmm.apply(linear_pose)
AddPyMOLObserver_to_conformation(linear_pose, True)

# 循环跑起来！
for i in range(100):
    trial_mover.apply(linear_pose)

# 获取能量最低的构象:
mc.recover_low(linear_pose)
mc.show_state()
mc.show_scores()

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.
protocols.moves.MonteCarlo: {0} MC: 1  24.7063  24.7063  24.7063  24.7063  0  0  0  accepted thermally
protocols.moves.TrialCounter: {0}            Shear trials=     46;  accepts= 0.6304;  energy_drop/trial=  -0.39129
protocols.moves.TrialCounter: {0}            Small trials=     54;  accepts= 0.7593;  energy_drop/trial=  -0.40653
protocols.moves.MonteCarlo: {0} MonteCarlo:: last_accepted_score,lowest_score: 24.7063 24.7063

结果好像并非我们所想象的那么顺利。最先开始的线性多肽并没有按照预期折叠成想要的形状。出现了什么问题？让我们看一下轨迹:

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

细心的同学可能会发现，我们模拟过程中没有处理侧链的结构，由于侧链没有进行变化，导致骨架发生变化时，侧链原子之间容易产生碰撞。导致每一步的能量基本都是升高的状态！这也是就在全原子打分函数下，能量面十分的崎岖，难以被遍历搜索。如果此时同时考虑侧链和主链，那自由度又过大，搜索效率下降。那在Rosetta中，这个问题是如何被解决的呢？那就是centroid模型与cen_std能量函数。

(图片来源: Leman J K, Weitzner B D, Lewis S M, et al. Macromolecular modeling and design in Rosetta: recent methods and frameworks[J]. Nature methods, 2020, 17(7): 665-680.)

二、粗粒化模型与能量面¶

在之前的章节中，我们提及过Rosetta中Pose有两种原子模型以及有粗粒化的打分函数。在粗粒化的原子模型和力场下，蛋白折叠的能量面平滑了许多，让我们对主链构象的搜索有了更大的移动空间。

(图片来源: Rosetta Workshop UNAM, by Possu Huang,Sergey Ovchinnikov)

在这一节，我们将尝试在粗粒化的打分函数的能量面上进行蛋白质折叠模拟。

In [5]:

# 转换pose、构建打分函数:
from pyrosetta.rosetta.protocols.simple_moves import SwitchResidueTypeSetMover

# 读入预先准备好的线性多肽;
linear_pose = pose_from_pdb('./data/linear_pep.pdb')

# 转换pose为粗粒化模型:
switch_mover = SwitchResidueTypeSetMover("centroid")
switch_mover.apply(linear_pose)

# cen_std打分函数
cen_score = create_score_function('cen_std')

# 重新定义MonteCarlo object:
mc = MonteCarlo(linear_pose, cen_score, kT)
trial_mover = TrialMover(rand_mover, mc)
pmm.apply(linear_pose)
AddPyMOLObserver_to_conformation(linear_pose, True)

core.import_pose.import_pose: {0} File './data/linear_pep.pdb' automatically determined to be of type PDB
core.chemical.GlobalResidueTypeSet: {0} Finished initializing centroid residue type set.  Created 69 residue types
core.chemical.GlobalResidueTypeSet: {0} Total time to initialize 0.022613 seconds.
basic.io.database: {0} Database file opened: scoring/score_functions/EnvPairPotential/env_log.txt
basic.io.database: {0} Database file opened: scoring/score_functions/EnvPairPotential/cbeta_den.txt
basic.io.database: {0} Database file opened: scoring/score_functions/EnvPairPotential/pair_log.txt
basic.io.database: {0} Database file opened: scoring/score_functions/EnvPairPotential/cenpack_log.txt

Out[5]:

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

In [6]:

# 循环跑起来！需要2-3min运行时间。
for i in range(5000):
    trial_mover.apply(linear_pose)

# 获取能量最低的构象:
mc.recover_low(linear_pose)
mc.show_state()
mc.show_scores()

# 可视化
linear_pose.dump_pdb('./data/centroid_mcmc.pdb')

protocols.moves.MonteCarlo: {0} MC: 1  9.81401  9.81401  9.81401  9.81401  0  0  0  accepted score beat last
protocols.moves.TrialCounter: {0}            Shear trials=   2501;  accepts= 0.9408;  energy_drop/trial=   0.00077
protocols.moves.TrialCounter: {0}            Small trials=   2499;  accepts= 0.8912;  energy_drop/trial=  -0.00199
protocols.moves.MonteCarlo: {0} MonteCarlo:: last_accepted_score,lowest_score: 9.81401 9.81401

Out[6]:

True

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

在构象搜索了5000步之后，我们发现在centroid原子模型和粗粒化的力场下，我们要模拟的多肽构象有了较大的构象变化，而不再是简单的线性结构。

思考: 为何还未成功折叠至我们的目标螺旋状态？

三、Fragment增强采样¶

最早在CASP8比赛中，Rosetta基于Fragment进行MCMC采样被成功地应用于蛋白结构预测(当然Fragment的方法做结构预测早已过时了)。这里使用的Fragment原理非常简单，首先通过序列的一级信息，可以预测得到相应的二级结构，只要从已有的17万个PDB结构数据中，找到类似序列且类似二级结构的多肽片段，将他们进行组装，就可以快速采样目标的结构。

Robetta的官方服务器上还留着Fragment生成的服务器，新手可以直接从这里获取Fragment。使用教育邮箱注册后就可以提交任务。此处最小的序列长度是27。

https://robetta.bakerlab.org/fragmentsubmit.jsp

在此处，我们更换一个更加复杂的蛋白作为案例: PDBID_1B72。

In [7]:

# 读取目标的pose
ref_pose = pose_from_pdb('./data/1b72.pdb')
ref_pose.sequence()  # 这就是1b72的序列，可以用这个去fragment服务器上提交任务。

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

Out[7]:

'LRTNFTTRQLTELEKEFHFNKYLSRARRVEIAATLELNETQVKIWFQNRRMKQKKRERE'

In [8]:

# 教案的制作中已经通过服务器准备好了相应的fragment文件，位于./data/*frags.
# 此处通过fragment相关的管理器mover读入fragment信息:
from pyrosetta.rosetta.core.fragment import ConstantLengthFragSet
from pyrosetta import pose_from_sequence
fragset3 = ConstantLengthFragSet(3)
fragset3.read_fragment_file("./data/3mer.frags")
fragset9 = ConstantLengthFragSet(9)
fragset9.read_fragment_file("./data/9mer.frags")

# 准备线性结构
frag_pose = pose_from_sequence('LRTNFTTRQLTELEKEFHFNKYLSRARRVEIAATLELNETQVKIWFQNRRMKQKKRERE', "centroid")

core.fragments.ConstantLengthFragSet: {0} finished reading top 200 3mer fragments from file ./data/3mer.frags
core.fragments.ConstantLengthFragSet: {0} finished reading top 200 9mer fragments from file ./data/9mer.frags

这里，我们不再使用SmallMover和ShearMover做主链结构的move，而是换用ClassicFragmentMover，这个mover可以将我们分离的3mer(3氨基酸片段)和9mer(9氨基酸片段)插入当前的线性结构当中。

In [9]:

#
from pyrosetta.rosetta.protocols.simple_moves import ClassicFragmentMover
from pyrosetta.rosetta.core.kinematics import MoveMap
movemap = MoveMap()
movemap.set_bb(True)

# 加载fragment插入mover
fragmover3 = ClassicFragmentMover(fragset3, movemap)
fragmover9 = ClassicFragmentMover(fragset9, movemap)

# 初始化combine mover
from pyrosetta.rosetta.protocols.moves import RandomMover
# RandomMover
rand_mover = RandomMover()
rand_mover.add_mover(fragmover3)
rand_mover.add_mover(fragmover9)

# 重新定义MonteCarlo object:
frag_mc = MonteCarlo(frag_pose, cen_score, kT)

from pyrosetta.rosetta.protocols.moves import TrialMover
trial_mover = TrialMover(rand_mover, frag_mc)

尝试运行以下用fragment插入法的mcmc代码。

In [10]:

# 循环跑起来！需要2-3min运行时间。
for i in range(10000):
    trial_mover.apply(frag_pose)

# 获取能量最低的构象:
frag_mc.recover_low(frag_pose)
frag_mc.show_state()
frag_mc.show_scores()

protocols.moves.MonteCarlo: {0} MC: 1  -15.9793  -15.9793  -15.9793  -15.9793  0  0  0  rejected
protocols.moves.TrialCounter: {0} ClassicFragmentM trials=  10000;  accepts= 0.0994;  energy_drop/trial=  -0.00288
protocols.moves.MonteCarlo: {0} MonteCarlo:: last_accepted_score,lowest_score: -15.9793 -15.9793

In [11]:

frag_pose.dump_pdb('./data/prediction.pdb')

Out[11]:

True

尝试比较预测结构和真实结构的RMSD差异:

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

思考1: 这种基于Fragment插入的结构预测方法有什么局限性？

思考2: 这种基于Fragment插入的结构预测方法最后采样得到的结构分布依然服从玻尔兹曼分布吗？为什么？

目前深度学习的方法预测蛋白质结构已经得到了长足的发展，已经全面超越了传统contacts+folding的建模策略。因此在结构预测领域已经很少使用。但目前Fragment的方法对于结构生成和蛋白质设计依然有着重要的作用。

练习题:¶

尝试构建一个每次在Fragment插入之后，接着使用SmallMover和ShearMover采样局部结构的结构搜索程序。【提示:可能需要2个mc object】

习题1中添加的模拟退火的步骤，查看对预测结构的影响？