分子对接(Moleculer-docking)理论

分子对接就是两个或多个分子之间通过几何匹配和能量匹配相互识别找到最佳匹配模式的过程。分子对接在酶学研究和药物设计中具有重要的应用意义。

分子对接计算是在受体活性位点区域通过空间结构互补和能量最小化原则来搜寻配体与受体是否能产生相互作用以及它们之间的最佳结合模式。其思想起源于Fisher的锁钥模型(即一把钥匙开一把锁),主要强调的是空间形状的匹配。但配体和受体的识别要比这个模型更加复杂。首先,配体和受体在对接过程中会由于相互适应而产生构象的变化。其次,分子对接还要求能量匹配,对接过程中结合自由能的变化决定了两个分子是否能够结合以及结合的强度。

1958年D.E.Koshland提出分子识别过程中的诱导契合概念,受体分子活性中心的结构原本并非与底物完全吻合,但其是柔软和可塑的。当配体与受体相遇时,可诱导受体构象发生相应的变化,从而便于他们的结合进而引起相应的反应。

分子对接方法根据不同的简化程度分为三类:刚性对接、半柔性对接和柔性对接。刚性对接指在对接过程中,受体和配体的构象不发生变化,适合研究比较大的体系如蛋白-蛋白之间以及蛋白-核酸之间,计算简单,主要考虑对象之间的契合程度。半柔性对接常用于小分子和大分子的对接,在对接过程中,小分子的构象可以在一定范围内变化,但大分子是刚性的。这样既可以在一定程度上考察柔性的影响,又能保持较高的计算效率。在药物设计和虚拟筛选过程中一般采用半柔性的分子对接方法。柔性对接方法一般用于精确研究分子之间的识别情况,由于允许对接体系的构象变化,可以提高对接准确性但耗时较长。

分子对接的目的是找到底物分子和受体分子最佳结合位置及其结合强度,最终可以获得配体和受体的结合构象,但这样的构象可以有很多,一般认为自由能最小的构象存在的概率最高。搜寻最佳构象就要用到构象搜索方法,常用的有系统搜索法和非系统搜索法。系统搜索法通过改变每个扭转角评估所有可能的结合构象,进而选取能量最低的。这一方法计算量非常大。因此通常使用非系统搜索法来寻找能量较低构象,常用方法有:分子动力学方法、随机搜索、遗传算法、距离几何算法等。随机搜索又包括完全随机算法、蒙特卡罗法和模拟退火法等。

分子对接软件包——AutoDock

AutoDock软件简介

AutoDOCK是一款由Scripps研究所的Olson实验室开发与维护的开源的分子模拟软件,最主要应用于蛋白与小分子配体的分子对接。其用户图形化界面(GUI)工具为AutoDOCK Tools(ADT)。AutoDock的新一代产品为AutoDock Vina。

AutoDock Vina使用拉马克遗传算法提高效率。软件把遗传算法和局部搜索结合在一起,遗传算法用于全局搜索,而局部搜索用于能量优化。为了加快计算速度,AutoDock Vina采用格点(grid)计算。首先在受体活性氨基酸附近划定一个长方体区域作为搜索空间,扫描不同类型的原子计算格点能量,在搜索空间内,调整配体的构象、位置和方向,进而评分、排序获得能量最低的构象作为输出结果。

对范德华相互作用的计算:每个格点上保存的范德华能量的值的数目与要对接的配体上的原子类型数目相同。如果一个配体中含有C、H、O三种原子类型,那么每个葛店需要用单个探针原子与来计算其与受体之间的范德华相互作用值。当配体与受体进行分子对接时,配体中某个原子和受体之间的相互作用能通过周围8个格点上的这种原子类型为探针的格点值用内插法得到。

静电相互作用的计算采用静电势格点。当配体与受体对接时,某个原子和受体之间的静电相互作用能通过周围格点上静电势以及原子上的部分电荷计算得到。

准备docking需要的受体(蛋白)和配体(化合物)

Docking算法需要每个原子带有电荷并且需要标记原子的属性。这些信息通常未包含在PDB文件中。我们需要在对蛋白和小分子的PDB文件预处理,生成同时包含以上信息和PDB文件中原子坐标信息的PDBQT格式文件。进一步地对于“柔性配体docking”,我们还需要定义配体的柔性部分和刚性部分。所有这些都可以通过软件AutoDock Tools (adt)来完成。

准备受体蛋白
  1. PDB文件(1hsg.pdb)中包含了蛋白、配体和水分子;首先提取出蛋白的坐标,即以关键字ATOM和TER开头的行存储到文件1hsg_prot.pdb。

    在Linux下,使用命令 egrep “^(ATOM|TER)” 1hsg.pdb > 1hsg_prot.pdb

  2. 启动AutoDockTools

    选择支持GUI的客户端,如MobaXterm,连接服务器。进入AutoDockTools安装目录的’bin’子目录,使用命令 ./adt &

  3. 依次点选 File–Read Molecule–1hsg_prot.pdb加载蛋白分子。

    ADT中按住左键拖动旋转分子结构;点击中键滚动缩放;按住右键移动晶体位置。

  4. 更改展示方式

    依次点选 Color–By Atom Type–All Geometries–OK。

  5. 加氢:晶体结构中通常缺少氢原子的坐标(因为氢原子电子少,且质子核对电子吸引能力弱,因此很难定位,具体见 http://www.uh.edu/~chembi/ChemSocRev_Jones_critical.pdf )。但是在docking过程中,氢原子尤其是极性氢原子对计算静电作用是必须的。因此我们需要给蛋白加上氢原子,之后氢原子会以白色短线形式出现。

    依次点选 Edit–Hydrogen–Add–Polar only–OK(只加极性H)

1hsg_prot_no_H.png
1hsg_prot_with_H.png

  1. 存储对蛋白的每个原子所做的修改和原子类型判断。

    依次点选 Grid–Macromolecule–Choose–1HSG_protein–Slect Molecule (ADT会弹出一个信息框包含程序所做的处理,比如合并非极性氢原子,计算原子局部电荷和判断原子类型,并提示保存Save–1hsg_prot.pdbqt)

  2. 在受体蛋白定义配体结合的3D搜索空间。如果我们事先不知道结合位点,理论上可以定义一个长方体盒子包含整个蛋白或者随便一个特定区域 (下文PDB文件解析中会提到PDB文件中有时会包含活性位点信息)。

    依次点选 Grid – Grid box ,将会在蛋白上画出一个长方体,并且有一个弹出框。在弹出框中,拖拽刻度线查看长方体的变化,完成设置。
    在这个例子中,我们知道结合位点,就选取以其为中心的一个小空间。
    设置Spacing (angstrom)为1埃 (这实际是一个换算系数, 相当于步长; 默认为0.375,是C-C单键长度的1/4,最大为1。spacing值与(各个维度上的点的数目+1)的乘积就是长方体Grid box的大小)。在我们调整的过程中,可以看到随着这个数值的变大,立方体也被放大了。
    另外我们设置x,y,z center为16,25,4,number of points in (x,y,z)-dimension为30,30,30(最大为126,必须为偶数,AutoDock会自动再每一维再加一个点)。记下我们设置的这些点,下面会用到。在刻度转盘处点击右键会弹出一个窗口,输入数字回车即可设置GRID的中心坐标和大小。较大的number of points in (xyz)-dimension和较小的Spacing会增加搜索的精度,同时需要花费更多的计算时间。

  3. 设置受体的柔性残基:

    在ADT中依次点选 Flexible Residues–Input–Choose Macromolecule–1hsg_prot; select–select from string–Residue: ARG8–Add–Dismiss, 8号ARG氨基酸残基就被选中了。
    再依次点选Flexible Residues–Choose Torsions in Currently Selected Residues将选择的残基标记为柔性残基并设置可扭转的数量。
    在分子显示窗口中分别点击两个残基上CA和CB原子之间的建,使之变为非扭转的(紫色显示)[Shift + 鼠标左键],这样两个残基中的32个键中有6个是可扭转的。这里设置配体的柔性残基或者使CA-CB的键为刚性都是可选操作。
    放在教程中只是用来展示怎么操作的,无其它指导意义。Flexible Residues–Output–Save Flexible PDBQT保存柔性残基文件。 Flexible Residues–Output–Save Rigid PDBQT保存柔性残基文件。

  4. 关掉grid和删除protein

    Grid Options–File–Close w/out saving; Edit–Delete–Delete Molecule–1hsg_prot–Continue

准备配体

与蛋白结构类似,配体的结构也缺少氢原子,我们需要添加氢原子并且定义哪些键是可以旋转的以用于柔性docking。

  1. 从PDB结构中提取配体的原子位置。indinavir的配体残基名字为MK1,以HETATM开头的行表示非核心多聚体的成分 (heteroatoms)(具体见PDB文件格式解释)。

    Linux系统下,运行 grep “^HETATM.*MK1” 1hsg.pdb > indinavir.pdb

  2. 将结构读入ADT:

    依次点选 File–Read Molecule–indinavir.pdb
    依次点选 Color–By Atom Type–All Geometreies–OK

  3. 给配体加上氢原子:

    Edit–Hydrogen–Add–Polar only–OK

  4. 在ADT中定义此化合物为配体,以便ADT为其计算局部电荷(partial charges)和设置可旋转配体键。

    依次点选 Ligand–input–Choose–indinavir–Select Molecule for AutoDock4
    这时会有一个弹出框显示ADT所做的操作,包括合并非极性氢(只在添加了的情况下)、计算电荷电量和设置旋转键。
    然后点选 Ligand–Output–Save as PDBQT 存储结果。

  5. 查看ADT检测出的旋转键:

    依次点选 Ligand–Torsion Tree–Choose Torsions,可以看到Number of rotatable bonds=14/32。

准备docking配置文件

docking配置文件包含了输入的受体(蛋白)、配体(化合物)和搜索参数的信息,为一个文本文件,名字任意,可以为conf.txt,内容如下:

receptor = 1hsg_prot.pdbqt
ligand = indinavir.pdbqt
num_modes = 50
out = dockingResult.pdbqt
log = docking.log
center_x = 16
center_y = 25
center_z = 4
size_x = 30
size_y = 30
size_z = 30
seed = 2009

receptor和ligand为输入文件的名字,与conf.txt在同一目录下; out为输出文件的名字;log为输出日志文件的名字。
centerhe和size定义搜索空间的位置和大小。
num_modes设置最多显示的结合模型,鉴于只输出符合能量值要求的结果,最后输出的结合模型数量可能少于这一数值。
seed设置随机数生成的种子,可以为任意整数。如果想重现之前的分析结果就需要使用相同的seed。

Docking 小分子化合物indinavir到HIV-1蛋白酶
  1. 使用AutoDock Vina执行docking预测

    在Linux终端执行 path_to_vina/vina –config conf.txt

  2. 输出结果包含两个文件,构象文件dockingResult.pdbqt和日志文件docking.log.

    dockingResult.pdbqt: 包含所有docking的模式,通常第一个为结合最好的构象,但如果前几个能量值相差不大时也有例外。
    docking.log: 日志文件,包含结合能量值(第一列,越低越稳定,默认由低到高排序,所以第一个为最好的构象)、每个构象与第一个构象的距离、每个构象与第一个构象的差别。

用Pymol可视化结果

打开PyMOL,依次点选File–Open文件类型选择All Files-选取结果dockingResult.pdbqt文件、原始蛋白和配体的pdb文件、原教程的pdbqt文件。

用Pymol展示配体和受体相互作用的原子和氢键

为了简化展示过程,我们设计了一个pml脚本 (脚本内有很详细的解释),只需要修改脚本里面受体和配体的名字,然后在PyMOL的命令行界面输入PyMOL> run display.pml即可获得展示结果。当然这个脚本也可以使用程序generatePmlForHbond.py生成。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
 ############################################################
###All one needs to do is replacing:                        ##
###  * Protein structure file: E:\docking\1hsg_prot.pdb     ##
###  * Protein name: 1hsg                                   ##
###  * Docking result file: E:\docking\indinavir.pdbqt      ##
###  * Docking result name (normally ligand name): indinavir##
 ############################################################
# The following 4 lines: 
	# 1. load protein structure and rename it
	# 2. add hydrogen (`h_add` uses a primitive algorithm to add hydrogens onto a molecule.)
	# 3. hide protein display
	# 4. show cartoon display for protein
load E:\yunpan\docking\1hsg_prot.pdb, 1hsg
h_add 1hsg
hide everything, 1hsg
show cartoon, 1hsg
cmd.spectrum("count", selection="1hsg", byres=1)

# The following 6 lines: 
	# 1. load ligand structure and rename it
	# 2. add hydrogen
	# 3. hide ligand display
	# 4. show ligand in sticks mode
	# 5. Set width of stick to 0.15
	# 6. Set atom color: C-white;N-blue;O-red
load E:\yunpan\docking\indinavir.pdbqt, indinavir
h_add indinavir
hide everything, indinavir
show sticks, indinavir
set stick_radius, 0.15
util.cbaw indinavir

# The following 1 line:
	# 1. Select metal ions
select metals, symbol mg+ca+fe+zn

# The following 2 lines:
	# 1. Set hydrogen donator
	# 2. Set hydrogen accrptor 
	# `select` creates a named selection from an atom selection. 
	# `select name, (selection)`
select h_donator,  (elem n,o and (neighbor hydro))
select h_acceptor, (elem o or (elem n and not (neighbor hydro)))

# The following 4 lines:
	# 1. Create link between ligand_h_acceptor and prot_h_donator  within given distance 3.2
	# 2. Create link between ligand_h_donator  and prot_h_acceptor within given distance 3.2
	#    Set filter 3.6 for ideal geometry and filter 3.2 for minimally acceptable geometry
	# 3. Set red color for ligand_h_acceptor and prot_h_donator 
	# 4. Set blue color for ligand_h_donator  and prot_h_acceptor
	# `distance` creates a new distance object between two selections. It will display all distances within the cutoff. Distance is also used to make hydrogen bonds like `distance hbonds, all, all, 3.2, mode=2`.
	# distance [ name [, selection1 [, selection2 [, cutoff [, mode ]]]]]
distance LaccPdon, (indinavir and h_acceptor), (1hsg and h_donator), 3.2
distance LdonPacc, (indinavir and h_donator), (1hsg and h_acceptor), 3.2
color red, LaccPdon
color blue, LdonPacc
#distance Fe_C20, (fep and name C20), (heme and name fe))

# The following 6 lines:
	# 1. Select non-hydro atoms of ligands
	# 2. Select protein atoms within 5A of selected atoms in last step
	# 3. Label alpha-c(ca) of selected residues with residue name and residue position
	# 4. Set label color back
	# 5. Set background white
	# 6. Hidden hydrogenes
select sele, indinavir & not hydro
select sele, byres (sele expand 5) & 1hsg
one_letter ={'VAL':'V', 'ILE':'I', 'LEU':'L', 'GLU':'E', 'GLN':'Q', \
'ASP':'D', 'ASN':'N', 'HIS':'H', 'TRP':'W', 'PHE':'F', 'TYR':'Y',    \
'ARG':'R', 'LYS':'K', 'SER':'S', 'THR':'T', 'MET':'M', 'ALA':'A',    \
'GLY':'G', 'PRO':'P', 'CYS':'C'}
label name ca & sele, "%s-%s" % (one_letter[resn],resi)
bg white
set label_color, black
hide (hydro)

# The follwing 5 lines
	# 1. Comment out this line
	# 2. Create an object `surrounding_res` to represent selected protein atoms
	#    `create`: creates a new molecule object from a selection. It can also be used to create states in an   existing object. 
	#    `create name, (selection)`
	# 3. Display created surface
	# 4. Set color for surrounding_res
	# 5. Set transparency for surrounding_res
	#    Transparency is used to adjust the transparency of Surfaces and Slices.    
	#    `set transparency, F, selection`
#show surface, 1hsg
create surrounding_res, sele
show surface, surrounding_res
color grey80, surrounding_res
set transparency, 0.5, surrounding_res

此外还可以使用如下脚本(list_hbonds.py)输出相互作用的原子及其位置。


# Copyright (c) 2010 Robert L. Campbell
from pymol import cmd

def list_hb(selection,selection2=None,cutoff=3.2,angle=55,mode=1,hb_list_name='hbonds'):
    """
    USAGE
  
    list_hb selection, [selection2 (default=None)], [cutoff (default=3.2)],
                       [angle (default=55)], [mode (default=1)],
                       [hb_list_name (default='hbonds')]
  
    The script automatically adds a requirement that atoms in the
    selection (and selection2 if used) must be either of the elements N or
    O.
  
    If mode is set to 0 instead of the default value 1, then no angle
    cutoff is used, otherwise the angle cutoff is used and defaults to 55
    degrees.
  
    e.g.
    To get a list of all H-bonds within chain A of an object
      list_hb 1abc & c. a &! r. hoh, cutoff=3.2, hb_list_name=abc-hbonds
  
    To get a list of H-bonds between chain B and everything else:
      list_hb 1tl9 & c. b, 1tl9 &! c. b
  
    """
    cutoff=float(cutoff)
    angle=float(angle)
    mode=float(mode)
    # ensure only N and O atoms are in the selection
    selection = selection + " & e. n+o"
    if not selection2:
        hb = cmd.find_pairs(selection,selection,mode=mode,cutoff=cutoff,angle=angle)
    else:
        selection2 = selection2 + " & e. n+o"
        hb = cmd.find_pairs(selection,selection2,mode=mode,cutoff=cutoff,angle=angle)
  
    # sort the list for easier reading
    hb.sort(lambda x,y:(cmp(x[0][1],y[0][1])))
  
    for pairs in hb:
        cmd.iterate("%s and index %s" % (pairs[0][0],pairs[0][1]), 'print "%1s/%3s`%s/%-4s " % (chain,resn,resi,name),')
        cmd.iterate("%s and index %s" % (pairs[1][0],pairs[1][1]), 'print "%1s/%3s`%s/%-4s " % (chain,resn,resi,name),')
        print "%.2f" % cmd.distance(hb_list_name,"%s and index %s" % (pairs[0][0],pairs[0][1]),"%s and index %s" % (pairs[1][0],pairs[1][1]))

#cmd.extend("list_hb",list_hb)
#if __name__ == "__main__":
cmd.load("E:/yunpan/docking/1hsg_prot.pdb", "1hsg")
cmd.h_add("(1hsg)")
cmd.load("E:/yunpan/docking/indinavir.pdbqt","indinavir")
cmd.h_add("(indinavir)")

h_donator  = "elem n,o & (neighbor hydro)"
h_acceptor = "elem o | (elem n & !(neighbor hydro))"

lacc = "indinavir & (elem o | (elem n & !(neighbor hydro)))"
ldon = "indinavir & (elem n,o & (neighbor hydro))"
pacc = "1hsg & (elem o | (elem n & !(neighbor hydro)))"
pdon = "1hsg & (elem n,o & (neighbor hydro))"

list_hb(ldon, pacc, hb_list_name="l2p_hbonds")
list_hb(lacc, pdon, hb_list_name="p2l_hbonds")

输出结果如下:
PyMOL>run E:/docking/list_hbonds.py
B/MK1’902/N4 B/GLY’27/O 3.03
B/MK1’902/O4 B/GLY’27/O 3.16
B/MK1’902/O2 A/ASP’25/OD1 2.77
B/MK1’902/O2 B/ASP’25/OD1 2.63

看上去比显示的氢键少了三个,这是因为我们在第二个函数中使用了H-键角度限制,如果在调用时给定参数list_hb(mode=0)则会获得一致结果。
after_run_list_hbond.png
after_run_hbond.pml.png
H-bond结果展示。第一张图为运行display.pml后的结果,蓝色虚线为氢键;第二张图为运行list_hbonds.py后的结果, 黄色虚线为氢键(覆盖了之前的蓝色)。可以通过点选LaccPdon, LdonPacc, l2p_hbonds显示不同的氢键。

展示疏水表面
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
# color_h
# -------
 
# PyMOL command to color protein molecules according to the Eisenberg hydrophobicity scale
 
#
# Source: http://us.expasy.org/tools/pscale/Hphob.Eisenberg.html
# Amino acid scale: Normalized consensus hydrophobicity scale
# Author(s): Eisenberg D., Schwarz E., Komarony M., Wall R.
# Reference: J. Mol. Biol. 179:125-142 (1984)
#
# Amino acid scale values:
#
# Ala:  0.620
# Arg: -2.530
# Asn: -0.780
# Asp: -0.900
# Cys:  0.290
# Gln: -0.850
# Glu: -0.740
# Gly:  0.480
# His: -0.400
# Ile:  1.380
# Leu:  1.060
# Lys: -1.500
# Met:  0.640
# Phe:  1.190
# Pro:  0.120
# Ser: -0.180
# Thr: -0.050
# Trp:  0.810
# Tyr:  0.260
# Val:  1.080
#
# Usage:
# color_h (selection)
#
from pymol import cmd
 
def color_h(selection='all'):
        s = str(selection)
	print s
        cmd.set_color('color_ile',[0.996,0.062,0.062])
        cmd.set_color('color_phe',[0.996,0.109,0.109])
        cmd.set_color('color_val',[0.992,0.156,0.156])
        cmd.set_color('color_leu',[0.992,0.207,0.207])
        cmd.set_color('color_trp',[0.992,0.254,0.254])
        cmd.set_color('color_met',[0.988,0.301,0.301])
        cmd.set_color('color_ala',[0.988,0.348,0.348])
        cmd.set_color('color_gly',[0.984,0.394,0.394])
        cmd.set_color('color_cys',[0.984,0.445,0.445])
        cmd.set_color('color_tyr',[0.984,0.492,0.492])
        cmd.set_color('color_pro',[0.980,0.539,0.539])
        cmd.set_color('color_thr',[0.980,0.586,0.586])
        cmd.set_color('color_ser',[0.980,0.637,0.637])
        cmd.set_color('color_his',[0.977,0.684,0.684])
        cmd.set_color('color_glu',[0.977,0.730,0.730])
        cmd.set_color('color_asn',[0.973,0.777,0.777])
        cmd.set_color('color_gln',[0.973,0.824,0.824])
        cmd.set_color('color_asp',[0.973,0.875,0.875])
        cmd.set_color('color_lys',[0.899,0.922,0.922])
        cmd.set_color('color_arg',[0.899,0.969,0.969])
        cmd.color("color_ile","("+s+" and resn ile)")
        cmd.color("color_phe","("+s+" and resn phe)")
        cmd.color("color_val","("+s+" and resn val)")
        cmd.color("color_leu","("+s+" and resn leu)")
        cmd.color("color_trp","("+s+" and resn trp)")
        cmd.color("color_met","("+s+" and resn met)")
        cmd.color("color_ala","("+s+" and resn ala)")
        cmd.color("color_gly","("+s+" and resn gly)")
        cmd.color("color_cys","("+s+" and resn cys)")
        cmd.color("color_tyr","("+s+" and resn tyr)")
        cmd.color("color_pro","("+s+" and resn pro)")
        cmd.color("color_thr","("+s+" and resn thr)")
        cmd.color("color_ser","("+s+" and resn ser)")
        cmd.color("color_his","("+s+" and resn his)")
        cmd.color("color_glu","("+s+" and resn glu)")
        cmd.color("color_asn","("+s+" and resn asn)")
        cmd.color("color_gln","("+s+" and resn gln)")
        cmd.color("color_asp","("+s+" and resn asp)")
        cmd.color("color_lys","("+s+" and resn lys)")
        cmd.color("color_arg","("+s+" and resn arg)")
cmd.extend('color_h',color_h)
 
def color_h2(selection='all'):
        s = str(selection)
	print s
	cmd.set_color("color_ile2",[0.938,1,0.938])
	cmd.set_color("color_phe2",[0.891,1,0.891])
	cmd.set_color("color_val2",[0.844,1,0.844])
	cmd.set_color("color_leu2",[0.793,1,0.793])
	cmd.set_color("color_trp2",[0.746,1,0.746])
	cmd.set_color("color_met2",[0.699,1,0.699])
	cmd.set_color("color_ala2",[0.652,1,0.652])
	cmd.set_color("color_gly2",[0.606,1,0.606])
	cmd.set_color("color_cys2",[0.555,1,0.555])
	cmd.set_color("color_tyr2",[0.508,1,0.508])
	cmd.set_color("color_pro2",[0.461,1,0.461])
	cmd.set_color("color_thr2",[0.414,1,0.414])
	cmd.set_color("color_ser2",[0.363,1,0.363])
	cmd.set_color("color_his2",[0.316,1,0.316])
	cmd.set_color("color_glu2",[0.27,1,0.27])
	cmd.set_color("color_asn2",[0.223,1,0.223])
	cmd.set_color("color_gln2",[0.176,1,0.176])
	cmd.set_color("color_asp2",[0.125,1,0.125])
	cmd.set_color("color_lys2",[0.078,1,0.078])
	cmd.set_color("color_arg2",[0.031,1,0.031])
        cmd.color("color_ile2","("+s+" and resn ile)")
        cmd.color("color_phe2","("+s+" and resn phe)")
        cmd.color("color_val2","("+s+" and resn val)")
        cmd.color("color_leu2","("+s+" and resn leu)")
        cmd.color("color_trp2","("+s+" and resn trp)")
        cmd.color("color_met2","("+s+" and resn met)")
        cmd.color("color_ala2","("+s+" and resn ala)")
        cmd.color("color_gly2","("+s+" and resn gly)")
        cmd.color("color_cys2","("+s+" and resn cys)")
        cmd.color("color_tyr2","("+s+" and resn tyr)")
        cmd.color("color_pro2","("+s+" and resn pro)")
        cmd.color("color_thr2","("+s+" and resn thr)")
        cmd.color("color_ser2","("+s+" and resn ser)")
        cmd.color("color_his2","("+s+" and resn his)")
        cmd.color("color_glu2","("+s+" and resn glu)")
        cmd.color("color_asn2","("+s+" and resn asn)")
        cmd.color("color_gln2","("+s+" and resn gln)")
        cmd.color("color_asp2","("+s+" and resn asp)")
        cmd.color("color_lys2","("+s+" and resn lys)")
        cmd.color("color_arg2","("+s+" and resn arg)")
cmd.extend('color_h2',color_h2)

将上面的脚本存储为color_h.py,在PyMOL界面运行File–Run–color_h.py,在命令行输入

PyMOl color_h –Show surface