在Pymol中获取所有二面角
我想在Pymol中获取一个蛋白质的所有二面角（包括phi、psi、chi1、chi2、chi3、chi4），但我只找到一个能显示phi和psi的函数。
比如说：
PyMOL>phi_psi 1a11
SER-2: ( 67.5, 172.8 )
GLU-3: ( -59.6, -19.4 )
LYS-4: ( -66.4, -61.7 )
MET-5: ( -64.1, -17.9 )
SER-6: ( -78.3, -33.7 )
THR-7: ( -84.0, -18.1 )
ALA-8: ( -85.7, -40.8 )
ILE-9: ( -75.1, -30.8 )
SER-10: ( -77.6, -47.0 )
VAL-11: ( -61.3, -27.4 )
LEU-12: ( -60.7, -47.5 )
LEU-13: ( -71.1, -38.6 )
ALA-14: ( -46.2, -50.7 )
GLN-15: ( -69.1, -47.4 )
ALA-16: ( -41.9, -52.6 )
VAL-17: ( -82.6, -23.7 )
PHE-18: ( -53.4, -63.4 )
LEU-19: ( -61.2, -30.4 )
LEU-20: ( -61.1, -32.3 )
LEU-21: ( -80.6, -60.1 )
THR-22: ( -45.9, -34.4 )
SER-23: ( -74.5, -47.8 )
GLN-24: ( -83.5, 11.0 )
但是缺少了手性角。有没有人知道怎么获取所有的二面角呢？
非常感谢！
4 个回答
这段内容来自于 计算Psi/Phi和/或二面角 Pymol - Pymol API
我一直找不到 /model/chain/segi/resn/resi/atom_name 这种原子选择方式。
直到我偶然发现了 选择宏：
选择宏可以让你用更简洁的方式表示一个很长的原子选择语句，比如：
PyMOL> select pept and segi lig and chain B and resi 142 and name CA
可以简化为：
PyMOL> select /pept/lig/B/142/CA
我试着计算 phy、psi 和 omega；仔细检查我的代码，看看我是否写对了：
#!/usr/bin/env python3
import pymol
from pymol import cmd, stored
pymol.finish_launching()
# pymol.finish_launching(['pymol', '-q', '-W' , '1200' , '-H' , '900' ])
cmd.fab('AAAGGAAGGVVPCAGCA', 'peptide')
# cmd.fab('AG', 'peptide')
# cmd.save('peptide.pdb' , 'peptide')
# cmd.hide('sticks' , 'peptide')
cmd.hide('everything' , 'peptide')
cmd.show('lines' , 'peptide')
a = cmd.select('selected' , 'peptide and resi 1-5 and name CA')
print('\n a ---> ' , a, '\n\n')
stored.dict = {}
cmd.iterate('selected' , "stored.dict[int(resi)] = [model, segi, chain, resn, int(resi)]")
print('\n stored.dict -------> ', stored.dict,'\n\n')
for i in stored.dict.keys() :
print( i , type(i))
print('stored.dict.keys ', stored.dict.keys() , type(stored.dict.keys()))
def dihedral_res( res_id , selected):
if int(res_id) - 1 in [i for i in selected.keys()] :
print(selected[res_id][3], ' can calculate phy', ' ____________' , res_id , int(res_id))
i = selected[res_id - 1]
# print(i)
s1 = "/{}/{}/{}/{}`{}/C".format(i[0] , i[1] , i [2] , i[3] , i[4])
# print(s1)
i = selected[res_id]
s2 = "/{}/{}/{}/{}`{}/N".format(i[0] , i[1] , i [2] , i[3] , i[4])
s3 = "/{}/{}/{}/{}`{}/CA".format(i[0] , i[1] , i [2] , i[3] , i[4])
s4 = "/{}/{}/{}/{}`{}/C".format(i[0] , i[1] , i [2] , i[3] , i[4])
cmd.select('dihedral_phy_'+str(res_id) , s1 +' '+ s2 +' '+ s3 +' '+ s4)
try:
dihedral_val_phy = cmd.get_dihedral(s1, s2, s3, s4, state=0)
except:
dihedral_val_phy = None
else :
print(selected[res_id][3], ' cannot calculate phy' , '____________' , res_id , int(res_id))
dihedral_val_phy = None
if int(res_id) + 1 in [i for i in selected.keys()] :
print(selected[res_id][3], ' can calculate psi', ' ____________' , res_id , int(res_id))
i = selected[res_id ]
# print(i)
s1 = "/{}/{}/{}/{}`{}/N".format(i[0] , i[1] , i [2] , i[3] , i[4])
# print(s1)
s2 = "/{}/{}/{}/{}`{}/CA".format(i[0] , i[1] , i [2] , i[3] , i[4])
s3 = "/{}/{}/{}/{}`{}/C".format(i[0] , i[1] , i [2] , i[3] , i[4])
i = selected[res_id + 1]
s4 = "/{}/{}/{}/{}`{}/N".format(i[0] , i[1] , i [2] , i[3] , i[4])
cmd.select('dihedral_psi_'+str(res_id) , s1 +' '+ s2 +' '+ s3 +' '+ s4)
try:
dihedral_val_psi = cmd.get_dihedral(s1, s2, s3, s4, state=0)
except:
dihedral_val_psi = None
else :
print(selected[res_id][3], ' cannot calculate psi' , '____________' , res_id , int(res_id))
dihedral_val_psi = None
if int(res_id) - 1 in [i for i in selected.keys()] :
print(selected[res_id][3], ' can calculate omega', ' ____________' , res_id , int(res_id))
i = selected[res_id - 1]
# print(i)
s1 = "/{}/{}/{}/{}`{}/CA".format(i[0] , i[1] , i [2] , i[3] , i[4])
s2 = "/{}/{}/{}/{}`{}/C".format(i[0] , i[1] , i [2] , i[3] , i[4])
# print(s1)
i = selected[res_id]
s3 = "/{}/{}/{}/{}`{}/N".format(i[0] , i[1] , i [2] , i[3] , i[4])
s4 = "/{}/{}/{}/{}`{}/CA".format(i[0] , i[1] , i [2] , i[3] , i[4])
cmd.select('dihedral_omega_'+str(res_id) , s1 +' '+ s2 +' '+ s3 +' '+ s4)
try:
dihedral_val_omega = cmd.get_dihedral(s1, s2, s3, s4, state=0)
except:
dihedral_val_omega = None
else :
print(selected[res_id][3], ' cannot calculate omega' , '____________' , res_id , int(res_id))
dihedral_val_omega = None
return dihedral_val_phy , dihedral_val_psi , dihedral_val_omega
for i in stored.dict.keys() :
print( dihedral_res( i, stored.dict) )
这是我计算一个假肽链 AAAGGAAGGVVPCAGCA 的前五个氨基酸的 phy、psi、omega 的代码输出的一部分。希望我方向没错。
ALA can calculate phy ____________ 3 3
ALA can calculate psi ____________ 3 3
ALA can calculate omega ____________ 3 3
(180.0, -179.39898681640625, -180.0)
这里是结果，突出显示了参与其中一个二面角的氨基酸：
对于其他的二面角，我想也可以用同样的方法来计算。我敢打赌有更好的方法来创建一个返回二面角的函数，但我对Python或蛋白质分子图形不太精通。为了完整起见，我附上了一张不同氨基酸侧链二面角命名规则的表格，希望这张表格有用 [取自 http://www.mlb.co.jp/linux/science/garlic/doc/commands/dihedrals.html ]：
当然可以！请看下面的内容：
在编程中，有时候我们需要让程序在特定的条件下执行某些操作。比如说，当用户点击一个按钮时，我们希望程序能做出反应。这种情况就需要用到“事件监听器”。简单来说，事件监听器就像一个守卫，它会一直关注某个特定的事件发生，比如鼠标点击、键盘输入等。
当事件发生时，事件监听器会立刻通知程序去执行相应的代码。这样，程序就能根据用户的操作做出反应，让整个应用变得更加互动和友好。
在编程中，我们通常会使用一些简单的语法来设置这些事件监听器。比如，我们可以用一行代码来告诉程序：“当这个按钮被点击时，请执行下面的代码。”这样，程序就会在用户点击按钮时，自动运行我们指定的操作。
总之，事件监听器是让程序能够响应用户操作的重要工具，它帮助我们创建更生动和有趣的应用体验。
#install bio3d library and call
library(bio3d)
#returns the file path of the current working directory.
getwd()
#sets the working directory to where you want.
setwd("home/R/Rscripts")
#fetches the pdb file from the protein data bank and saves to dataframe 'pb'
pb <- read.pdb("insert PDB ID")
#trim to protein only
pb.prot <- trim.pdb(pb, "protein")
#calculates the torsion angles of the protein and save to dataframe 'tor'
tor <- torsion.pdb(pb.prot)
#to get the number of rows and columns of 'tor'
dim(tor$tbl)
#identify each row by their chain, residue ID and residue Number obtained from your PDB entry
res_label <- paste(pb.prot$atom$chain[pb.prot$calpha], pb.prot$atom$resid[pb.prot$calpha], pb.prot$atom$resno[pb.prot$calpha], sep="-")
rownames(tor$tbl) <- res_label
#creates a table of the torsion angle
torsion <- tor$tbl
#For example, to look at the angles for VAL, residue 223 from chain A
tor$tbl["A-GLY-65",]
#convert "double" into a datatype
dataframe_df=as.data.frame.matrix(torsion)
#write dataframe to a .csv file
write.csv(dataframe_df, file="name.csv", row.names=TRUE,col.names=TRUE)
你可以很简单地在R语言中做到这一点。这里有一个链接，里面有关于如何计算主链和侧链扭转角的信息：http://thegrantlab.org/bio3d/html/torsion.pdb.html
但首先你需要安装Bio3D这个R语言的包：http://thegrantlab.org/bio3d/download
安装好这个包后，在R的控制台输入library(bio3d)来加载它。
>library(bio3d)
这段R脚本可以帮助你解决问题：
#返回当前工作目录的文件路径。
getwd()
#设置工作目录到你想要的位置。
setwd("home/R/Rscripts")
#从蛋白质数据银行获取pdb文件，并保存到数据框'pb'中。
pb <- read.pdb("insert PDB ID")
#只保留蛋白质部分。
pb.prot <- trim.pdb(pb, "protein")
#计算蛋白质的扭转角，并保存到数据框'tor'中。
tor <- torsion.pdb(pb.prot)
#获取'tor'的行数和列数。
dim(tor$tbl)
#通过链、残基ID和残基编号来识别每一行，这些信息来自你的PDB条目。
res_label <- paste(pb.prot$atom$chain[pb.prot$calpha], pb.prot$atom$resid[pb.prot$calpha], pb.prot$atom$resno[pb.prot$calpha], sep="-")
rownames(tor$tbl) <- res_label
#创建一个扭转角的表格。
torsion <- tor$tbl
#例如，查看来自链A的VAL，残基223的角度。
tor$tbl["A-VAL-223",]
#将表格写入一个文件。
write.table(torsion, file = "torsion_angles.txt", quote = F, sep = "\t")
你在工作目录中保存的输出文件将包含链、残基ID、残基编号及其对应的phi、psi、chi1、chi2、chi3、chi4和chi5值的表格。祝你好运！
你可以通过 get_dihedral 获取任意的二面角。首先，创建四个选择，每个选择里放一个原子，然后像这样使用：
get_dihedral s1, s2, s3, s4
这个功能在 Python 接口中叫做 cmd.get_dihedral()。我建议你写一个 Python 脚本，利用这个函数和 cmd.iterate() 来遍历每个残基。你可以创建一个字典，这样在处理每个残基时，就能查找出一组原子四元组，这些四元组用来定义卡角（chi-angles）。