Rosetta的能量函数
Rosetta能量函数的详细情况在我的另外一篇文章里，这里就不过多的赘述了，下面将具体到pyrosetta中去了解能量函数的使用。
pyrosetta中能量函数的使用
#导入pyrosetta
import pyrosetta
pyrosetta.init()
from pyrosetta.rosetta.core.scoring import *
#导入pdb的一条链
ras = pyrosetta.pose_from_pdb("6Q21_A.pdb")
#导入能量函数
from pyrosetta.teaching import *
sfxn = get_score_function(True) # ref2015 all-atom energy function
为了得到一个蛋白的能量得分，首先要通过pyrosetta.teaching 里面导入get_score_function(is_fullatom: bool) 函数，当为True时将返回ref2015 all-atom energy function，当为False时将返回centroid score function。在Centroid模型下，氨基酸的侧链被一个有一定半径大小的粗粒化球所替代，而在full-atom模型下，所有的侧链原子被显式的表示出来。
同时，我们这也可以通过改变不同能量项的权重去改变能量函数。
sfxn2 = ScoreFunction()
sfxn2.set_weight(fa_atr, 1.0)
sfxn2.set_weight(fa_rep, 1.0)
# 这里我们的能量函数只包括范德华力，并将范德华力的引力（fa_atr）和斥力（fa_rep）的权重都设置成1。
然后我们可以同时运行上面定义的两个能量函数sfxn和sfxn2来对pose进行操作。
print(sfxn(ras))
print(sfxn2(ras))
basic.io.database: Database file opened: scoring/score_functions/elec_cp_reps.dat
core.scoring.elec.util: Read 40 countpair representative atoms
core.pack.dunbrack.RotamerLibrary: shapovalov_lib_fixes_enable option is true.
core.pack.dunbrack.RotamerLibrary: shapovalov_lib::shap_dun10_smooth_level of 1( aka lowest_smooth ) got activated.
core.pack.dunbrack.RotamerLibrary: Binary rotamer library selected: /Users/kathyle/Computational Protein Prediction and Design/PyRosetta4.Release.python36.mac.release-208/pyrosetta/database/rotamer/shapovalov/StpDwn_0-0-0/Dunbrack10.lib.bin
core.pack.dunbrack.RotamerLibrary: Using Dunbrack library binary file '/Users/kathyle/Computational Protein Prediction and Design/PyRosetta4.Release.python36.mac.release-208/pyrosetta/database/rotamer/shapovalov/StpDwn_0-0-0/Dunbrack10.lib.bin'.
core.pack.dunbrack.RotamerLibrary: Dunbrack 2010 library took 0.450999 seconds to load from binary
1215.729069796814
154.59159174026854
sfxn.show(ras) #显示出能量函数的每一项
core.scoring:
------------------------------------------------------------
Scores Weight Raw Score Wghtd.Score
------------------------------------------------------------
fa_atr 1.000 -1039.246 -1039.246
fa_rep 0.550 1193.837 656.611
fa_sol 1.000 682.582 682.582
fa_intra_rep 0.005 700.419 3.502
fa_intra_sol_xover4 1.000 46.564 46.564
lk_ball_wtd 1.000 -14.597 -14.597
fa_elec 1.000 -195.387 -195.387
pro_close 1.250 97.210 121.513
hbond_sr_bb 1.000 -41.656 -41.656
hbond_lr_bb 1.000 -28.352 -28.352
hbond_bb_sc 1.000 -13.111 -13.111
hbond_sc 1.000 -7.771 -7.771
dslf_fa13 1.250 0.000 0.000
omega 0.400 41.525 16.610
fa_dun 0.700 1296.642 907.650
p_aa_pp 0.600 -25.496 -15.298
yhh_planarity 0.625 0.000 0.000
ref 1.000 47.114 47.114
rama_prepro 0.450 197.781 89.002
---------------------------------------------------
Total weighted score: 1215.729
如果我们想要计算每一个氨基酸的能量，
ras.energies().show() #返回pose中每一个氨基酸的能量以及各项的能量
ras.energies().show(24) #返回24号氨基酸的能量以及各项的能量
根据之前的文章我们知道，范德华力的吸引力和排斥力、溶剂化能以及静电力都是与原子的类型和就离相关的，所以我们可以根据etable_atom_pair_energies(res1, atom_index_1, res2, atom_index_2, sfxn)两个原子来计算出他们之间的力。
res24 = ras.residue(24)
res20 = ras.residue(20)
res24_atomN = res24.atom_index("N")
res20_atomO = res20.atom_index("O")
pyrosetta.etable_atom_pair_energies(res24, res24_atomN, res20, res20_atomO, sfxn)
(-0.1505855046001568, 0.0, 0.5903452111877215, 2.173111777247698)
pose.energies().total_energies()[fa_atr] # 返回所有能量中fa_atr的贡献
pose.energies().residue_total_energies(5)[fa_atr] # 返回第五个氨基酸残基能量中fa_atr的贡献
# 如果出现NameError: name 'fa_atr' is not defined，记得导入from pyrosetta.rosetta.core.scoring import *
hbond_set = hbonds.HBondSet() # 创建一个set来储存氢键的能量和信息
pose.update_residue_neighbors() # 基于相邻的残基更新pose的energies
hbonds.fill_hbond_set(pose, False,hbond_set)
hbond_set = pyrosetta.rosetta.core.pose.Pose.get_hbonds(pose) # 将hbonds填充到hbond_set 里面
hbond_set.show(pose) # 展示出pose中所有的氢键以及能量
emap = EMapVector() # 创建一个energy map
res102 = pose.pdb_info().pdb2pose("D", 102)
res408 = pose.pdb_info().pdb2pose("A", 408)
# 获得两个氨基酸在pose中的number信息
sfxn.eval_ci_2b(pose.residue(res102), pose.residue(res408), pose, emap) # 计算这两个氨基酸之间的能量并将其储存在energy map里面
emap[fa_atr] #返回 energymap 中的fa_atr
