问题描述
我们创建了一个函数 superimpose，它接受 7 个参数（PDB_id_1、Chain_id_1、Res_1、PDB_id_2、Chain_id_2、Res_2、输出文件名）。
然后，我们使用 get_coordinates 函数从两个文件的输入中提到的残基中提取了各个原子的各自坐标。
之后，我们使用了 RMSD 解并将其添加到 superimpose 函数中以计算 RMSD 旋转和平移矩阵。
然后我们将函数 download_pdb 合并到 superimpose 函数中，以便下载所需的 PDB 文件。
最后，我们通过函数 apply_rot_trans
脚本是：
import Bio
from Bio.PDB import *
from numpy import *
from numpy.linalg import svd,det
from pathlib import Path
def superimpose(PDB_id_1,Chain_id_1,Res_1,PDB_id_2,Chain_id_2,Res_2,out_filename):
# Parse PDB file
inputs = []
coordinates_both_files = []
inputs.append([PDB_id_1,Res_1])
inputs.append([PDB_id_2,Res_2])
for i in range(len(inputs)):
# downloading the pdb file
filename = download_pdb(inputs[i][0])
p=PDBParser()
s=p.get_structure(inputs[i][0],filename)
# Getting the required Chain
# hierarchy: structure[model][chain][residue no.][atom]
chain = s[0][inputs[i][1]]
coordinates_both_files.append(get_coordinates(chain,inputs[i][2]))
##################### superimpose ######################
# Nr of atoms
x = coordinates_both_files[0]
y = coordinates_both_files[1]
N=x.shape[1]
########################
# Center x and y
x=center(x)
y=center(y)
# correlation matrix
#r=y*x.T
r = x*y.T
# SVD of correlation matrix
#v,s,wt=svd(r)
v,wt=svd(r)
#w=wt.T
w=wt.T
#vt=v.T
vt=v.T
# Rotation matrix
u=w*vt
# Check for roto-reflection
if det(u)<0:
z=diag((1,1,-1))
u=w*z*vt
s[2]=-s[2]
# Translation matrix
#t = y - (u * x)
t = u * x+y
print("Translation Matrix:",t)
# Calculate RMSD
#e0=sum(x.A*x.A+y.A*y.A)
e0=sum(y.A*y.A+x.A*x.A)
print("E0 ",e0)
rmsd=sqrt((e0-2*sum(s))/N)
print('RMSD (svd) ',rmsd)
print("u estimated ")
print(u)
# Calculate RMSD from the coordinates
d=x-u*y
d=d.A*d.A
rmsd=sqrt(sum(d)/N)
print('RMSD (real) ',rmsd)
############ applying rotation and translation to all atoms of 1BBH ############
apply_rot_trans(u,t,inputs[0][0],out_filename)
def get_coordinates(chain,res_id_list):
backbone = ["N","CA","C","O"]
coordinates = []
for res_id in res_id_list:
for atom in backbone:
try:
# hierarchy: structure[model][chain][residue no.][atom]
# since we have a chain we can get the atom: atom = chain[residue no.][atom]
atom=chain[res_id][atom]
c=atom.get_coord()
coordinates.append(c)
except:
pass
# Turn coordinate list in 3xn numpy matrix
coordinates=matrix(coordinates) # nx3
coordinates=coordinates.T # 3xn
return coordinates
def center(m):
# Returns centered m
# Calculate center of mass of x
center_of_mass_m=m.sum(1)/m.shape[1]
# Center m
centered_m=m-center_of_mass_m
return centered_m
def download_pdb(pdb_id):
out_dir = Path.cwd()
filename = "pdb"+pdb_id+".ent"
# try and except block to prevent re downloading the same file on multiple executions of the code
try:
# trying to open the file
# will print the given error message if PDB file already exists
fr = open(filename,"r")
fr.close()
print("pdb file already exist")
except:
# will download the file only if it doen not already exists
pdbl = PDBList()
pdbl.retrieve_pdb_file(pdb_id,pdir = out_dir,file_format = "pdb")
return filename
def apply_rot_trans(u,pdb_id,out_filename):
superimposed_coord = []
filename = "pdb"+pdb_id+".ent"
p=PDBParser()
s=p.get_structure(pdb_id,filename)
atom_list = list(s.get_atoms())
for a in atom_list:
a.coord = dot(u[1],a.coord)+t
print(a.coord)
# writing the output to a pdb file
io = PDBIO()
struct = io.set_structure(s)
io.save(struct)
PDB_id_1 = "1BBH"
Chain_id_1 = "B"
Res_1 = [16,121,124,125]
PDB_id_2 = "5GYR"
Chain_id_2 = "B"
Res_2 = [16,125]
out_filename = "rotation.pdb"
superimpose(PDB_id_1,out_filename)
但它给了我以下错误：
TypeError Traceback (most recent call last)
<ipython-input-16-e0ae1f929075> in <module>()
150 Res_2 = [16,125]
151 out_filename = "rotation.pdb"
--> 152 superimpose(PDB_id_1,out_filename)
<ipython-input-16-e0ae1f929075> in superimpose(PDB_id_1,out_filename)
79
80 ############ applying rotation and translation to all atoms of 1BBH ############
---> 81 apply_rot_trans(u,out_filename)
82
83 def get_coordinates(chain,res_id_list):
<ipython-input-16-e0ae1f929075> in apply_rot_trans(u,out_filename)
140 io = PDBIO()
141 struct = io.set_structure(s)
--> 142 io.save(struct)
143
144
~\Anaconda3\lib\site-packages\Bio\PDB\PDBIO.py in save(self,file,select,write_end,preserve_atom_numbering)
376 resseq,377 icode,--> 378 chain_id,379 )
380 fp.write(s)
~\Anaconda3\lib\site-packages\Bio\PDB\PDBIO.py in _get_atom_line(self,atom,hetfield,segid,atom_number,resname,resseq,icode,chain_id,charge)
227 charge,228 )
--> 229 return _ATOM_FORMAT_STRING % args
230
231 else:
TypeError: only size-1 arrays can be converted to Python scalars
解决方法
暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！
如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@）