Notebook

Atom & Residue Layer¶

@Author: 吴炜坤

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

在Rosetta中原子和残基的处理方式和众多的分子模拟工具相类似。基于拓扑结构定义以及参数实例化的方式高效地处理每个“独立构象的”氨基酸残基。这些独立的氨基酸残基串起来，就构成了整个蛋白质的构象。

一、AtomType和ResidueType¶

AtomType和ResidueType其实可以理解为模板文件，它定义了一个原子的基本信息以及氨基酸残基的基本拓扑信息。(定义残基应该长什么样子)

AtomType: 记录了原子的元素名、电荷、LJ半径，LK溶剂体积等一些最基本的参数，这些参数可以在Rosetta的源代码中被找到。
ResidueType: 记录了残基的三字母缩写，单字母缩写，组成残基的原子类型(AtomType的信息，如原子名，原子电荷)，原子间的成键和配对方式(Bond connectivity)，原子间的内坐标二面角信息（internal coordinates）, 这些信息都被储存在Rosetta params文件里。

此处以THR残基的Params文件做范例:

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

In [1]:

# 初始化PyRosetta
from pyrosetta import *
init()

PyRosetta-4 2021 [Rosetta PyRosetta4.conda.mac.cxx11thread.serialization.python37.Release 2021.31+release.c7009b3115c22daa9efe2805d9d1ebba08426a54 2021-08-07T10:04:12] 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 r292 2021.31+release.c7009b3115c c7009b3115c22daa9efe2805d9d1ebba08426a54 http://www.pyrosetta.org 2021-08-07T10:04:12
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=1698516215 seed_offset=0 real_seed=1698516215 thread_index=0
basic.random.init_random_generator: {0} RandomGenerator:init: Normal mode, seed=1698516215 RG_type=mt19937

二、Pose Residue¶

2.1 从ResidueType生成Residue¶

有了氨基酸的模板信息，我们只要对氨基酸模板中的二面角参数进行“赋值”，就可以用这个氨基酸来描述PDB文件中任意相同类型氨基酸的具体构象了。而其他的信息如化学键的连接顺序，原子类型，原子电荷都是固有属性，无法被变更。

Rosetta处理的逻辑:

根据Pose中的残基名，找到对应的ResidueType(模板)，复制这个模板; 根据PDB中的坐标信息，计算得到每个4原子间的二面角，并将这些信息填入ICOOR_INTERNAL(内坐标地)信息当中; 根据新的ResidueType信息，生成一个独立的Residue对象。

更具体地:

Rosetta中根据ResidueTypeSets来储存这些模板，我们第一步，需要先把对应的ResidueType(模板找出来);
通过ResidueFactory来制造和生产一个标准的Residue对象
将Residue加载到Pose中。

In [2]:

# 通过ResidueFactory创建一个新的Residue(ALA)对象:
from pyrosetta.rosetta.core.chemical import ChemicalManager
from pyrosetta.rosetta.core.conformation import ResidueFactory

chm = ChemicalManager.get_instance()
residue_type_sets = chm.residue_type_set("fa_standard")

# 使用残基名来匹配residue_type模板:
ala_residue_type = residue_type_sets.name_map('ALA')

# 通过ResidueFactory制造一个Residue:
new_residue= ResidueFactory.create_residue(residue_type_sets.name_map('ALA')) # 创建一个residue object

print(new_residue)

core.chemical.GlobalResidueTypeSet: {0} Finished initializing fa_standard residue type set.  Created 983 residue types
core.chemical.GlobalResidueTypeSet: {0} Total time to initialize 0.649089 seconds.
Residue 0: ALA (ALA, A):
Base: ALA
 Properties: POLYMER PROTEIN CANONICAL_AA ALIPHATIC METALBINDING ALPHA_AA L_AA
 Variant types:
 Main-chain atoms:  N    CA   C  
 Backbone atoms:    N    CA   C    O    H    HA 
 Side-chain atoms:  CB  1HB  2HB  3HB 
Atom Coordinates:
   N  : 0, 0, 0
   CA : 1.458, 0, 0
   C  : 2.00885, 1.42017, 0
   O  : 1.38304, 2.33899, -0.528697
   CB : 1.9878, -0.772763, -1.19909
   H  : -0.491969, 0.763905, -0.44104
   HA : 1.79666, -0.489627, 0.91315
  1HB : 3.07772, -0.76375, -1.18511
  2HB : 1.63286, -1.80244, -1.15412
  3HB : 1.63327, -0.30698, -2.11716
Mirrored relative to coordinates in ResidueType: FALSE

In [3]:

# 创建一个1个氨基酸的pose:
from pyrosetta import pose_from_sequence
pose = pose_from_sequence('A')
print(pose)

PDB file name: A
Total residues: 1
Sequence: A
Fold tree:
FOLD_TREE  EDGE 1 1 -1

In [4]:

# 在第一个氨基酸后添加一个ALA，并使用理想构象
pose.append_polymer_residue_after_seqpos(new_residue, seqpos=1, build_ideal_geometry=True) 
print(pose)

PDB file name: A
Total residues: 2
Sequence: AA
Fold tree:
FOLD_TREE  EDGE 1 1 -1  EDGE 1 2 -1

这个过程中从ResidueFactory生成的Residue其实都是理想型的构象，因此这过程中还涉及到调整 $\chi$ ， $\phi$ ， $\psi$ ， $\omega$ 角的过程，这一个内容在下一节进行介绍。

2.2 索引Pose获取Residue对象¶

如果现在已经生成好了一个PDB相应的Pose了，我们可以通过索引的机制去将其中的每一个Reside对象进行获取。

In [5]:

# 初始化
from pyrosetta import init, pose_from_pdb
init()
pose = pose_from_pdb('./data/4R80.clean.pdb')

PyRosetta-4 2021 [Rosetta PyRosetta4.conda.mac.cxx11thread.serialization.python37.Release 2021.31+release.c7009b3115c22daa9efe2805d9d1ebba08426a54 2021-08-07T10:04:12] 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 r292 2021.31+release.c7009b3115c c7009b3115c22daa9efe2805d9d1ebba08426a54 http://www.pyrosetta.org 2021-08-07T10:04:12
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=-156941142 seed_offset=0 real_seed=-156941142 thread_index=0
basic.random.init_random_generator: {0} RandomGenerator:init: Normal mode, seed=-156941142 RG_type=mt19937
core.import_pose.import_pose: {0} File './data/4R80.clean.pdb' automatically determined to be of type PDB
core.conformation.Conformation: {0} [ WARNING ] missing heavyatom:  OXT on residue SER:CtermProteinFull 76
core.conformation.Conformation: {0} [ WARNING ] missing heavyatom:  OXT on residue SER:CtermProteinFull 152

In [6]:

# 通过Pose和residue的Pose编号，我们可以直接获取到Residue，并提取里面的具体信息:
residue1 = pose.residue(1) 
residue2 = pose.residue(2) 
print(residue1)

Residue 1: PRO:NtermProteinFull (PRO, P):
Base: PRO
 Properties: POLYMER PROTEIN CANONICAL_AA LOWER_TERMINUS ALIPHATIC METALBINDING ALPHA_AA L_AA
 Variant types: LOWER_TERMINUS_VARIANT
 Main-chain atoms:  N    CA   C  
 Backbone atoms:    N    CA   C    O    HA 
 Side-chain atoms:  CB   CG   CD   NV  CAV  1HB  2HB  1HG  2HG  1HD  2HD  1H   2H  
Atom Coordinates:
   N  : 35.432, -0.708, 7.647
   CA : 35.959, 0.478, 8.332
   C  : 36.62, 1.469, 7.374
   O  : 36.946, 1.11, 6.24
   CB : 36.987, -0.1, 9.317
   CG : 37.119, -1.563, 8.97
   CD : 35.846, -1.957, 8.305
   NV : 35.4257, -0.707336, 7.64284 (virtual)
  CAV : 35.8996, 0.441347, 8.25002 (virtual)
   HA : 35.141, 0.97715, 8.8721
  1HB : 37.9438, 0.433685, 9.21849
  2HB : 36.6417, 0.0476102, 10.3509
  1HG : 37.9855, -1.72036, 8.31089
  2HG : 37.3003, -2.15428, 9.87969
  1HD : 36.0441, -2.75894, 7.57861
  2HD : 35.1232, -2.2907, 9.06406
  1H  : 35.7595, -0.700397, 6.70024
  2H  : 34.4274, -0.649029, 7.64284
Mirrored relative to coordinates in ResidueType: FALSE

通过Residue这个对象，我们可以很方便地获取到每个氨基酸的坐标残基以及部分原子级别的信息。

2.3 Residue Name¶

在Rosetta中，处于N端和C端的氨基酸通常都有质子化的现象(NH3+或COO-)，因此和处于肽链中的氨基酸残基的拓扑结构不同(原子数目也不同), 如果需要描述这些质子化的氨基酸信息，就要在原有的params文件基础上添加新的原子连接信息。如果每一种氨基酸都有2种质子化状态，就至少需要3x20个params文件！

为了高效描述和减少Params的文件数量，Rosetta提出了Patch的机制(补丁)，通过补丁系统，对于同一种氨基酸的不同状态只需要在原有的params文件内容中把额外的内容临时加上即可（类比打补丁）。

如形成了质子化的N端氨基酸是以"被修饰"的状态进行处理, 因此他们的残基命名也带上了补丁字样，比如，1号氨基酸名称为PRO:NtermProteinFull，其中NtermProteinFull就是他的"补丁名"。补丁的内容记录在NtermProteinFull这个patch文件中。

因此Rosetta中的残基名信息有好几种形式存在:

annotated_name: 带有Patch信息的残基名
name1: 单字母缩写
name3: 三字母缩写

In [7]:

# 获取残基的Rosetta残基名、单字母缩写、三字母缩写、标注名
annotated_name = pose.residue(1).annotated_name()
name1 = residue1.name1()
name3 = residue1.name3()
print(annotated_name)
print(name1)
print(name3)

P[PRO:NtermProteinFull]
P
PRO

2.4 Residue Property¶

每种氨基酸残基都有自带的一些性质/属性标签。这些字段由Params文件中的property字段记录。通过residue对象的Properties输出的行，也可以轻松获取这些信息。

In [8]:

# 注意Properties的输出信息: POLYMER PROTEIN CANONICAL_AA LOWER_TERMINUS ALIPHATIC METALBINDING ALPHA_AA L_AA
print(residue1)

Residue 1: PRO:NtermProteinFull (PRO, P):
Base: PRO
 Properties: POLYMER PROTEIN CANONICAL_AA LOWER_TERMINUS ALIPHATIC METALBINDING ALPHA_AA L_AA
 Variant types: LOWER_TERMINUS_VARIANT
 Main-chain atoms:  N    CA   C  
 Backbone atoms:    N    CA   C    O    HA 
 Side-chain atoms:  CB   CG   CD   NV  CAV  1HB  2HB  1HG  2HG  1HD  2HD  1H   2H  
Atom Coordinates:
   N  : 35.432, -0.708, 7.647
   CA : 35.959, 0.478, 8.332
   C  : 36.62, 1.469, 7.374
   O  : 36.946, 1.11, 6.24
   CB : 36.987, -0.1, 9.317
   CG : 37.119, -1.563, 8.97
   CD : 35.846, -1.957, 8.305
   NV : 35.4257, -0.707336, 7.64284 (virtual)
  CAV : 35.8996, 0.441347, 8.25002 (virtual)
   HA : 35.141, 0.97715, 8.8721
  1HB : 37.9438, 0.433685, 9.21849
  2HB : 36.6417, 0.0476102, 10.3509
  1HG : 37.9855, -1.72036, 8.31089
  2HG : 37.3003, -2.15428, 9.87969
  1HD : 36.0441, -2.75894, 7.57861
  2HD : 35.1232, -2.2907, 9.06406
  1H  : 35.7595, -0.700397, 6.70024
  2H  : 34.4274, -0.649029, 7.64284
Mirrored relative to coordinates in ResidueType: FALSE

In [9]:

# 通过Residue Object还可以判断氨基酸化学性质,残基归类属性
print(residue1.is_polar())
print(residue1.is_aromatic())
print(residue1.is_charged())
print(residue1.is_DNA())
print(residue1.is_protein())
print(residue1.is_aramid())

False
False
False
False
True
False

运行结果表明，Residue1是非极性非带电的蛋白质氨基酸。

2.5 原子信息索引¶

氨基酸的原子id信息也是由Params文件所定义，因此有固定的顺序, 并且从1开始编号。因此获取原子的信息，第一件需要做的事是索引原子号:

通过原子id进行索引
通过原子名进行索引

In [10]:

# 通过原子名获取:
ca_id = residue1.atom_index('CA') # 残基内atom id
print(ca_id)

In [11]:

# 获取总原子数量:
natoms = residue1.natoms()
print(f'total atoms: {natoms}')

# 通过原子id获取原子信息:
for atom_id in range(1, natoms+1):
    print(atom_id, residue1.atom_name(atom_id))

total atoms: 18
1  N  
2  CA 
3  C  
4  O  
5  CB 
6  CG 
7  CD 
8  NV 
9 CAV 
10  HA 
11 1HB 
12 2HB 
13 1HG 
14 2HG 
15 1HD 
16 2HD 
17 1H  
18 2H

2.6 原子坐标¶

获取原子坐标的方式有两种:

通过residue的xyz函数获取（封装，更方便）
通过原子对象获取

In [12]:

# 直接通过原子名获取原子的xyz坐标。
ca_xyz = pose.residue(1).xyz("CA")
print(ca_xyz)

# 直接通过原子号获取原子的xyz坐标。
atom1_xyz = pose.residue(1).xyz(1)
print(atom1_xyz)

      35.95900000000000      0.4780000000000000       8.332000000000001
      35.43200000000000     -0.7080000000000000       7.647000000000000

有了原子的id信息，就可以获取到具体的某个原子对象:

通过原子名获取对象
通过原子id获取对象

In [13]:

# 先获取原子对象，再获取坐标。
# 1 原子名获取对象
atom_ca = residue1.atom("CA")

# 2 通过原子id获取对象
atom1 = residue1.atom(1)

通过原子对象，也可以获取坐标信息:

In [14]:

print(atom_ca.xyz())
print(atom1.xyz())

      35.95900000000000      0.4780000000000000       8.332000000000001
      35.43200000000000     -0.7080000000000000       7.647000000000000

2.7 原子类型信息的获取¶

每个原子的信息由AtomType模板进行定义，类似于ResidueType, 我们可以获取每个原子的细节信息，如原子的Rosetta类型、原子的元素名、范德华半径、Lazaridis Karplus溶剂化参数等.

通过atom_type方法，就可以迅速获取原子的性质参数。

In [15]:

# 获取1号原子的atomtype所有的信息:
atom_type1 = residue1.atom_type(1)
print(atom_type1)

Atom Type: Nlys
	element: N
	Lennard Jones: radius=1.80245 wdepth=0.161725
	Lazaridis Karplus: lambda=3.5 volume=16.514 dgfree=-20.8646
	properties: DONOR 
Extra Parameters: 1.75 1.55 0.79 1.55 1.44 1.5 1.55 -20 -10.695 -1.145 -20 -0.62 0 0 0 1.85 8.52379 0.025 0.01 0.005 -289.292 -0.697267 -1933.88 -1.56243 -93.2613 93.2593 0.00202205 715.165 74.6559 -74.6539 0.00268963 -1282.36 0.633 -0.367 0.926 -0.537 0.633 -0.367

In [16]:

# 获取Rosetta的atomtype_name
atomtype_name = atom_type1.atom_type_name()
print(atomtype_name)

Nlys

In [17]:

# 元素名
atom_element = atom_type1.element()
print(atom_element)

In [18]:

# LJ参数
lj_radius = atom_type1.lj_radius()
lk_volume = atom_type1.lk_volume()
lk_dgfree = atom_type1.lk_dgfree()

print(lj_radius, lk_volume, lk_dgfree)

1.802452 16.514 -20.864641

思考¶

如果现在要定一个全新的非天然氨基酸，我都需要哪些具体的信息？

In [ ]: