Error in Saving Files

Question

[anasjamshed]

I have made this script:

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)):
        # 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, 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 

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,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

[Ram]

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.

[anasjamshed]

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

[Ram]

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.

[anasjamshed]

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

[Ram]

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.