PyMol - script

从大二开始接触PyMol至今已有4年了，之前只会用PyMol展示PDB文件。记得有一回需要找蛋白质所有催化残基3.6埃以内的活性位点。想到PyMol打开PDB文件时肯定已经完成了解析，就想写一个脚本拿到感兴趣的位点然后遍历所有原子坐标找到3.6埃以内的残疾。由于对PyMol不是很熟悉，这个想法在脑海中闪现之后就夭折了。后来我写了个C++程序，从PDB文件开始解析，真是废了九牛二虎之力，现在想想，当时真是太simple，还是要学习啊。

Pymol调用用户脚本

今天学习了一下PyMol调用脚本的过程，于是有了这篇博客。

PyMol是用Python实现的开源软件。用户可以用Python写一个自定义函数（没试过类，以后有需要再试试），关键是如何让让自己写的函数嵌入到PyMol。

主要过程如下（流程图）：

Created with Raphaël 2.1.0run scriptcmd.extendEnd

所有操作在PyMol命令窗口中完成。

这里新建了两个文件：pymolrc.pml，start.py。两个文件都放在宿主目录下。请不要修改第一个文件的名字，第二个文件的名字请自便。如果修改的话，下面也要相应修改成你的文件名。Windows的宿主目录位于C:/Users/ste/，C代表系统盘的位置。Linux宿主目录不知道的话，敲个cd直接回车，进入的目录就是你的宿主目录。PyMol在启动时会自动调用宿主目录下的pymolrc文件，在该文件中添加

run start.py

假设在start.py中定义了一个名为hello的函数，并且不传入任何参数。接下来使用

cmd.extend("hello", hello)

将hello函数注册为PyMol命令。当然你也可以将上面这行代码添加到start.py中，接下来就可以像使用label、create这些内置PyMol命令一样使用hello了

command形式：

hello

PyMol API形式：

hello()

使用上述两个文件有一个好处：所有自定义命令会在PyMol启动时注册为PyMol命令，从而扩展了PyMol命令及API。

下面分别是我的pymolrc.pml和start.py文件

set label_position, (3,2,1)
run ~/start.py

cd E:\Dropbox\endnote\CRISPR.Data\PDB

#http://pymolwiki.org/index.php/Simple_Scripting#Python_Version_Supported_by_PyMOL
oneLetter ={'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', 'ADE':'A', 'URI':'U', 'GUA':'G', 'CYT':'C', 'THY':'T'}
mylabel = 'oneLetter[resn] + resi'

RPOTEIN = 1
NUCLEICACID = 2
def mylabel(sel, mode=1):
mode = int(mode)
if mode == 1:
cmd.label(sel, 'oneLetter[resn] + resi')
else:
cmd.label(sel, 'resn + resi')

cmd.extend("mylabel", mylabel)

from pymol import stored

def getseq2(sel, mode = RPOTEIN):
stored.seq = []
mode = int(mode)
# iterate over state 1, or the sel -- this just means
# for each item in the selection do what the next parameter says.
# And, that is to append the residue name to the stored.seq
# array.
if mode == RPOTEIN:
sel= sel + " and n. CA"
cmd.iterate_state(1, selector.process(sel), "stored.seq.append(oneLetter[resn])")
else:
sel1= sel + " and n. p"
sel2= sel + " and n. c1'"
print cmd.count_atoms(sel1)
sel = sel1 if cmd.count_atoms(sel1) > cmd.count_atoms(sel2) else sel2
cmd.iterate_state(1, selector.process(sel), "stored.seq.append(resn)")
print sel
print ''.join(stored.seq)

cmd.extend("getseq2", getseq2)

#iterate al atoms in selection, extract resn per resi
def getseq(sel):
stored.seq = []
atoms = cmd.get_model(sel)
#get first atom
atom = atoms.atom[0]
resiCurr = atom.resi
if atom.resn in oneLetter.keys():
stored.seq.append(oneLetter[atom.resn])
elif len(atom.resn) == 1:
stored.seq.append(atom.resn)
else:
stored.seq.append(atom.resn + ' ')

for atom in atoms.atom:
if atom.resi != resiCurr:
resiCurr = atom.resi
if atom.resn in oneLetter.keys():
stored.seq.append(oneLetter[atom.resn])
elif len(atom.resn) == 1:
stored.seq.append(atom.resn)
else:
stored.seq.append(' ' + atom.resn + ' ')
print re.sub('(.{60})', r'\1\n', ''.join(stored.seq))

cmd.extend("getseq", getseq)

def mytest():
atoms = cmd.get_model("chain A and resi 2-3 and n. C+O+N+CA")
for at in atoms.atom:
print "ATOM DEFINITION: "+at.chain+" "+at.resn+" "+str(at.resi)+" " +at.name+" "+str(at.index)+" "+str(at.b)+" "+str(at.coord[0])+" "+str(at.coord[1])+" "+str(at.coord[2])

cmd.extend("mytest", mytest)

def myexit():
os._exit(0)
cmd.extend("myexit", myexit)

def getseqbychain():
chains = cmd.get_chains('all')
for chain in chains:
#print chain
sel = 'chain %s' % chain
print '>%s' % chain
getseq(sel)
cmd.extend('getseqbychain', getseqbychain)

操作原子

到此为止，用户就可以顺利向PyMol注册自定义命令了。接下来要做的就是如何从PyMol获得原子坐标等基本信息。

核心思路是通过selection获得一个object，这个object包含原子信息。

可以参考这篇wiki：get_model

我的start.py有两个函数可以参考：getseq，getseqbychain。getseq主要是根据传入的selection获得所有atom instance，然后遍历之。因为一个残基含有多个原子，这里需要过滤。同一个残基上的原子含有相同的resi，根据这一点去除冗余。

刚开始想到的是使用残基的一个特征原子来代替残基，从而避免遍历所有的原子。对蛋白质，可以使用CA；对核酸就比较复杂了。有些核酸的某个残基没有C1’或者P，导致selection漏掉了这个残基，函数getseq2就是按照这种思路实现的。

资源推荐

pymolwiki

这里可以找到你想要的PyMol命令的用法及示例。

simple scripting

pymol mail list

PyMol用户众多，你遇到的问题，别人很可能已经遇到，而且很可能已经有了解决方案。

Google

Google在任何时候都是一个很好的助手！