Error in Saving Files
0
0
Entering edit mode
18 months ago

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,Chain_id_1, Res_1])
inputs.append([PDB_id_2,Chain_id_2, Res_2])
for i in range(len(inputs)):
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, s, 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

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()

except:
pdbl = PDBList()
pdbl.retrieve_pdb_file(pdb_id, pdir = out_dir, file_format = "pdb")

return filename

def apply_rot_trans(u,t,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,121,124,125]
out_filename = "rotation.pdb"
superimpose(PDB_id_1, Chain_id_1, Res_1, PDB_id_2, Chain_id_2, Res_2, out_filename)


But is giving me following error :

TypeError                                 Traceback (most recent call last)
<ipython-input-16-e0ae1f929075> in <module>()
150 Res_2 = [16,121,124,125]
151 out_filename = "rotation.pdb"
--> 152 superimpose(PDB_id_1, Chain_id_1, Res_1, PDB_id_2, Chain_id_2, Res_2, out_filename)

<ipython-input-16-e0ae1f929075> in superimpose(PDB_id_1, Chain_id_1, Res_1, PDB_id_2, Chain_id_2, Res_2, out_filename)
79
80     ############ applying rotation and translation to all atoms of 1BBH ############
---> 81     apply_rot_trans(u,t,inputs[0][0],out_filename)
82
83 def get_coordinates(chain, res_id_list):

<ipython-input-16-e0ae1f929075> in apply_rot_trans(u, t, pdb_id, 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


Due to this error, I am unable to save the structure file. Plz help me

python PDB • 513 views
2
Entering edit mode

You see how you're trying different things and running into different errors? That's how you learn. Keep at it, you'll figure out the root cause by yourself soon. Talk to a programmer close to you if you'd like to expedite the process. An online forum is not the best place for "debug my code" style questions. Even StackOverflow would not allow these types of questions, and we are not a programming forum.

By the way, this is the third question you've started for the exact same topic. This is not good etiquette as you should keep your discussion to one thread. I don't want to discourage you, but this question might be closed as a duplicate, and opening new questions for the same discussion might get you suspended, so please be more mindful of your behavior.

0
Entering edit mode

the topic is not the same. I am finding errors in saving file

0
Entering edit mode

Are you telling me it's a different problem because a different line went wrong in your code? Unless you had working code that broke when you made a change and after fixing it, broke in a different way when you made a different change, it does not count as separate topics. You're writing code from scratch that's running into problems, that's one topic - not two independent bugs.

I'd recommend using base python instead of ipython notebook, and even switching to python 3 (you seem to be using python2). The two factors above can only add to the complication, and you will need to debug as close to the interpreter as you can.

0
Entering edit mode

It's python 3, not Python 2. I am used to with Jupiter notebook and thonny

0
Entering edit mode

If you can account for factors that the IDE and Jupyter introduce, debug with them. Debugging run time errors is always simpler without an IDE, in my opinion. I did not realize Python 3 supported return xyz statements, sorry for assuming you're using python 2.