Pyrosetta入门

---

前言

最近又使用Rosetta的需求，故今天想写一下Pyrosetta的入门学习笔记。

---

一、安装

conda config --add channels https://levinthal:paradox@conda.graylab.jhu.edu
conda install pyrosetta

二、蛋白相关处理

1.处理PDB

# 读入一个pdb结构文件
from pyrosetta import *
pose = pose_from_pdb("inputs/5tj3.pdb")
pose.sequence()
pose.annotated_sequence()


# 处理（清理）一个pdb文件
from pyrosetta.toolbox import cleanATOM
cleanATOM("inputs/5tj3.pdb")
pose_clean = pose_from_pdb("inputs/5tj3.clean.pdb")
pose_clean.sequence()
pose_clean.annotated_sequence()

2.残基信息

# 获得残基的名字
print(pose.total_residue())
residue20 =  pose.residue(20)
print(residue20.name())

# 判断是不是金属
zn_resid = pose.pdb_info().pdb2pose('A', 601)
res_zn = pose.residue(zn_resid) 
res_zn.is_metal() 

3.获取空间特征

这里的角度还是没太搞懂。。待补充，也许是二面角吗

# 获得残基空间中的角度的方法
resid = pose.pdb_info().pdb2pose('A', 28)
print("phi:", pose.phi(resid))
print("psi:", pose.psi(resid))
print("chi1:", pose.chi(1, resid))

# 获得键长的两种方法
# 方法一：
res_28 = pose.residue(resid)
N28 = AtomID(res_28.atom_index("N"), resid)
CA28 = AtomID(res_28.atom_index("CA"), resid)
C28 = AtomID(res_28.atom_index("C"), resid)
print(pose.conformation().bond_length(N28, CA28))
print(pose.conformation().bond_length(CA28, C28))

# 方法二：
N_xyz = res_28.xyz("N")
CA_xyz = res_28.xyz("CA")
C_xyz = res_28.xyz("C")
N_CA_vector = CA_xyz - N_xyz
CA_C_vector = CA_xyz - C_xyz
print(N_CA_vector.norm())
print(CA_C_vector.norm())

# 获得键角的方法
angle = pose.conformation().bond_angle(N28, CA28, C28)
print(angle)

tripeptide = pose_from_sequence("AAA")

orig_phi = tripeptide.phi(2)
orig_psi = tripeptide.psi(2)
print("original phi:", orig_phi)
print("original psi:", orig_psi)
print("xyz coordinates:", tripeptide.residue(2).xyz("CB"))

tripeptide.set_phi(2, -60)
tripeptide.set_psi(2, -43)
print("new phi:", tripeptide.phi(2))
print("new psi:", tripeptide.psi(2))
print("xyz coordinates:", tripeptide.residue(2).xyz("CB"))

三、能量函数

1.分数函数基础

import pyrosetta
ras = pyrosetta.pose_from_pdb("inputs/6Q21_A.pdb")

from pyrosetta.teaching import *
sfxn = get_score_function(True)
print(sfxn)

这样会看到项，权重和能量方法选择。

2.自定义能量函数

sfxn2 = ScoreFunction()
sfxn2.set_weight(fa_atr, 1.0)
sfxn2.set_weight(fa_rep, 1.0)

print(sfxn(ras))
print(sfxn2(ras))

可以通过设置能量项前的权重自定义能量函数。

3.能量分解

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

将能量分解到残基和原子：

# 将能量分解到残基
print(ras.energies().show(24))

# 将能量分解到原子对（输出的是两原子之间的吸引、排斥、溶剂化和静电）
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)

分解两个残基之间的相互作用的能量：

from pyrosetta import *
sfxn = get_score_function(True)

pose = pyrosetta.toolbox.pose_from_rcsb("1YY9")
res102 = pose.pdb_info().pdb2pose("D", 102)
res408 = pose.pdb_info().pdb2pose("A", 408)

emap = EMapVector()
sfxn.eval_ci_2b(pose.residue(res102), pose.residue(res408), pose, emap)
print(emap[fa_atr])
print(emap[fa_rep])
print(emap[fa_sol])

---

总结

未完待续。