Merge pull request #31 from Pranavkhade/hdanm_improvement

Hdanm improvement

packman/__init__.py

@@ -15,4 +15,4 @@
 
 
 #VERSION CHANGE HERE CHANGES IT IN docs AND setup.py
-__version__='1.4.3'
+__version__='1.4.4'

packman/anm/hd_anm.py

@@ -40,6 +40,7 @@
     * Pranav Khade(https://github.com/Pranavkhade)
 """
 
+import logging
 from .. import molecule, Atom, Model, Protein
 from ..constants import amino_acid_molecular_weight
 from ..constants import atomic_weight
@@ -551,7 +552,7 @@ def calculate_cross_correlation( self, n_modes = "all" ):
 
                 self.crosscorrelation_matrix[i][j] = trace_H_inv_ij / numpy.sqrt( trace_H_inv_ii*trace_H_inv_jj )
 
-    def calculate_movie(self, mode_number, scale=1.5, n=20, extrapolation="curvilinear", ftype='cif'):
+    def calculate_movie(self, mode_number, scale=1.5, n=20, extrapolation="curvilinear", ftype='cif', ca_to_aa=False):
         """This function generates the dynamic 3D projection of the normal modes obtained using hd-ANM. The 3D projection can be linearly extrapolated or curvilinearly extrapolated depending on the choices. The first frame is the original structure and the projection progresses in positive (+) direction, returning to original structure and then in negative direction (-) again returning to the original structure.
 
         Args:
@@ -560,6 +561,7 @@ def calculate_movie(self, mode_number, scale=1.5, n=20, extrapolation="curviline
             n (int)                             : Number of frames in output (should be =>8 and ideally multiple of 4)   Defaults to 20
             extrapolation (linear/curvilinear)  : Extrapolation method                                                   Defaults to "curvilinear"
             ftype (string)                      : Extension of the output file (.cif / .pdb)
+            ca_to_aa (boolean)                  : If only single atom type is used (often C-alpha atom only), enabling extrapolates the C-alpha motion to all the atoms. (Default: False)
 
         Note:
             - Scale and n parameters should be redesigned.
@@ -569,6 +571,11 @@ def calculate_movie(self, mode_number, scale=1.5, n=20, extrapolation="curviline
             True if successful; false otherwise.
         """
 
+        if(ca_to_aa):
+            if(len(list(set([i.get_name() for i in self.atoms]))) > 1 ): logging.warning('Multiple atom types were detected, but the "ca_to_aa" option is still "True." If the model is not coarse-grained, consider making it "False." See the function description for more details.')
+        else:
+          if(len(list(set([i.get_name() for i in self.atoms]))) == 1 ): logging.info('A single type of atom is detected. Try enabling the "ca_to_aa" parameter. See the function description for more details.')
+
         x0=numpy.array([i.get_location() for i in self.atoms])
         d0=[i.get_domain_id() for i in self.atoms]
         new_coords=[]
@@ -579,19 +586,39 @@ def calculate_movie(self, mode_number, scale=1.5, n=20, extrapolation="curviline
             for j in movement:
                 HingeResidue=0
                 AtomsOfTheFrame = {}
+                non_model_atoms = 0 #Only when ca_to_aa option is enabled.
                 for numi,i in enumerate(x0):
                     if(d0[numi][0]=='D'):
                         D_delta_phi= self.eigen_vectors[:,mode_number][self.domain_info[d0[numi]][0]*6 : (self.domain_info[d0[numi]][0]*6)+6]
                         D_COM=self.domain_info[d0[numi]][1]
                         new_x=i[0]+scale*j* ( D_delta_phi[0] + (D_delta_phi[4]*(i[2]-D_COM[2])) - (D_delta_phi[5]*(i[1]-D_COM[1])) )
                         new_y=i[1]+scale*j* ( D_delta_phi[1] + (D_delta_phi[3]*(i[2]-D_COM[2])) - (D_delta_phi[5]*(i[0]-D_COM[0])) )
                         new_z=i[2]+scale*j* ( D_delta_phi[2] + (D_delta_phi[3]*(i[1]-D_COM[1])) - (D_delta_phi[4]*(i[0]-D_COM[0])) )
+                        if(ca_to_aa):
+                            for x in self.atoms[numi].get_parent().get_atoms():
+                                temp_new_x=x.get_location()[0]+scale*j* ( D_delta_phi[0] + (D_delta_phi[4]*(x.get_location()[2]-D_COM[2])) - (D_delta_phi[5]*(x.get_location()[1]-D_COM[1])) )
+                                temp_new_y=x.get_location()[1]+scale*j* ( D_delta_phi[1] + (D_delta_phi[3]*(x.get_location()[2]-D_COM[2])) - (D_delta_phi[5]*(x.get_location()[0]-D_COM[0])) )
+                                temp_new_z=x.get_location()[2]+scale*j* ( D_delta_phi[2] + (D_delta_phi[3]*(x.get_location()[1]-D_COM[1])) - (D_delta_phi[4]*(x.get_location()[0]-D_COM[0])) )
+
+                                non_model_atoms+=1
+                                AtomsOfTheFrame[len(x0)+non_model_atoms] = Atom(x.get_id() , x.get_name(), numpy.array([temp_new_x,temp_new_y,temp_new_z]), x.get_occupancy(), x.get_bfactor(), x.get_element(), x.get_charge(), x.get_parent() )
+
                     if(d0[numi][0]=='H'):
                         delta=self.eigen_vectors[:,mode_number][6*len(self.domain_info):][HingeResidue*3:(HingeResidue*3)+3]
                         new_x=i[0]+scale*j*delta[0]
                         new_y=i[1]+scale*j*delta[1]
                         new_z=i[2]+scale*j*delta[2]
                         HingeResidue+=1
+                        if(ca_to_aa):
+                            for x in self.atoms[numi].get_parent().get_atoms():
+                                temp_new_x=x.get_location()[0]+scale*j*delta[0]
+                                temp_new_y=x.get_location()[1]+scale*j*delta[1]
+                                temp_new_z=x.get_location()[2]+scale*j*delta[2]
+
+                                non_model_atoms+=1
+                                AtomsOfTheFrame[len(x0)+non_model_atoms] = Atom(x.get_id() , x.get_name(), numpy.array([temp_new_x,temp_new_y,temp_new_z]), x.get_occupancy(), x.get_bfactor(), x.get_element(), x.get_charge(), x.get_parent() )
+
+                    #All the atoms in the models (CA or all atom); Atoms from ca_to_aa option are not included.
                     new_x , new_y, new_z = new_x.real , new_y.real , new_z.real
                     new_coords.append([new_x,new_y,new_z])
                     currentatm = Atom(self.atoms[numi].get_id() , self.atoms[numi].get_name(), numpy.array([new_x,new_y,new_z]), self.atoms[numi].get_occupancy(), self.atoms[numi].get_bfactor(), self.atoms[numi].get_element(),self.atoms[numi].get_charge(), self.atoms[numi].get_parent() )
@@ -611,6 +638,7 @@ def calculate_movie(self, mode_number, scale=1.5, n=20, extrapolation="curviline
             for j in movement:
                 HingeResidue=0
                 AtomsOfTheFrame = {}
+                non_model_atoms = 0 #Only when ca_to_aa option is enabled.
                 for numi,i in enumerate(x0):
                     if(d0[numi][0]=='D'):
                         D_delta_phi= self.eigen_vectors[:,mode_number][self.domain_info[d0[numi]][0]*6 : (self.domain_info[d0[numi]][0]*6)+6]
@@ -633,14 +661,44 @@ def calculate_movie(self, mode_number, scale=1.5, n=20, extrapolation="curviline
                         + ( D_mu[0]*D_mu[2]*(1-numpy.cos(scale*j*Q_D_n)) - D_mu[1]*numpy.sin(scale*j*Q_D_n) ) * (i[0] -D_COM[0]) \
                         + ( D_mu[1]*D_mu[2]*(1-numpy.cos(scale*j*Q_D_n)) + D_mu[0]*numpy.sin(scale*j*Q_D_n) ) * (i[1] -D_COM[1]) \
                         + ( numpy.cos(scale*j*Q_D_n)+ D_mu[2]**2 * (1-numpy.cos(scale*j*Q_D_n)) )             * (i[2] -D_COM[2])
-
+
+                        if(ca_to_aa):
+                            for x in self.atoms[numi].get_parent().get_atoms():
+                                temp_new_x= D_COM[0]+ scale*j*D_delta_phi[0] \
+                                + ( numpy.cos(scale*j*Q_D_n)+ D_mu[0]**2 * (1-numpy.cos(scale*j*Q_D_n)) )             * (x.get_location()[0] -D_COM[0]) \
+                                + ( D_mu[0]*D_mu[1]*(1-numpy.cos(scale*j*Q_D_n)) - D_mu[2]*numpy.sin(scale*j*Q_D_n) ) * (x.get_location()[1] -D_COM[1]) \
+                                + ( D_mu[0]*D_mu[2]*(1-numpy.cos(scale*j*Q_D_n)) + D_mu[1]*numpy.sin(scale*j*Q_D_n) ) * (x.get_location()[2] -D_COM[2])
+
+                                temp_new_y= D_COM[1]+ scale*j*D_delta_phi[1] \
+                                + ( D_mu[0]*D_mu[1]*(1-numpy.cos(scale*j*Q_D_n)) + D_mu[2]*numpy.sin(scale*j*Q_D_n) ) * (x.get_location()[0] -D_COM[0]) \
+                                + ( numpy.cos(scale*j*Q_D_n)+D_mu[1]**2 * (1-numpy.cos(scale*j*Q_D_n)) )              * (x.get_location()[1] -D_COM[1]) \
+                                + ( D_mu[1]*D_mu[2]*(1-numpy.cos(scale*j*Q_D_n)) - D_mu[0]*numpy.sin(scale*j*Q_D_n) ) * (x.get_location()[2] -D_COM[2])
+
+                                temp_new_z= D_COM[2]+ scale*j*D_delta_phi[2] \
+                                + ( D_mu[0]*D_mu[2]*(1-numpy.cos(scale*j*Q_D_n)) - D_mu[1]*numpy.sin(scale*j*Q_D_n) ) * (x.get_location()[0] -D_COM[0]) \
+                                + ( D_mu[1]*D_mu[2]*(1-numpy.cos(scale*j*Q_D_n)) + D_mu[0]*numpy.sin(scale*j*Q_D_n) ) * (x.get_location()[1] -D_COM[1]) \
+                                + ( numpy.cos(scale*j*Q_D_n)+ D_mu[2]**2 * (1-numpy.cos(scale*j*Q_D_n)) )             * (x.get_location()[2] -D_COM[2])
+
+                                non_model_atoms+=1
+                                AtomsOfTheFrame[len(x0)+non_model_atoms] = Atom(x.get_id() , x.get_name(), numpy.array([temp_new_x,temp_new_y,temp_new_z]), x.get_occupancy(), x.get_bfactor(), x.get_element(), x.get_charge(), x.get_parent() )
+
 
                     if(d0[numi][0]=='H'):
                         delta=self.eigen_vectors[:,mode_number][6*len(self.domain_info):][HingeResidue*3:(HingeResidue*3)+3]
                         new_x=i[0]+scale*j*delta[0]
                         new_y=i[1]+scale*j*delta[1]
                         new_z=i[2]+scale*j*delta[2]
                         HingeResidue+=1
+
+                        if(ca_to_aa):
+                            for x in self.atoms[numi].get_parent().get_atoms():
+                                temp_new_x=x.get_location()[0]+scale*j*delta[0]
+                                temp_new_y=x.get_location()[1]+scale*j*delta[1]
+                                temp_new_z=x.get_location()[2]+scale*j*delta[2]
+
+                                non_model_atoms+=1
+                                AtomsOfTheFrame[len(x0)+non_model_atoms] = Atom(x.get_id() , x.get_name(), numpy.array([temp_new_x,temp_new_y,temp_new_z]), x.get_occupancy(), x.get_bfactor(), x.get_element(), x.get_charge(), x.get_parent() )
+
                     new_x , new_y, new_z = new_x.real , new_y.real , new_z.real
                     new_coords.append([new_x,new_y,new_z])
                     currentatm = Atom(self.atoms[numi].get_id() , self.atoms[numi].get_name(), numpy.array([new_x,new_y,new_z]), self.atoms[numi].get_occupancy(), self.atoms[numi].get_bfactor(), self.atoms[numi].get_element(),self.atoms[numi].get_charge(), self.atoms[numi].get_parent() )

packman/apps/hdanm.py

@@ -26,7 +26,7 @@ def generate_output():
         numpy.savetxt("eigenvectors.csv", Model.get_eigenvectors(), delimiter=",")
 
         for i in range(6,6+args.modes,1):
-            Model.calculate_movie(i,scale=args.scale,n=args.frames)
+            Model.calculate_movie(i,scale=args.scale,n=args.frames,ca_to_aa=args.ca_to_aa)
 
     if args.make_tar:
 
@@ -44,4 +44,4 @@ def generate_output():
                 tar.close()
 
     else:
-        generate_output()
+        generate_output()

packman/bin/GUI.py

@@ -14,6 +14,7 @@
 #                                             GUI                                                #
 ##################################################################################################
 '''
+from xmlrpc.client import boolean
 import numpy
 import logging
 
@@ -27,10 +28,28 @@
     from tkinter.messagebox import showinfo, showerror
     from tkinter import N, S, E, W, Grid, END, BooleanVar, StringVar, DISABLED
     from tkinter.filedialog import askopenfilename
+    import codecs
+    import os.path
 except:
     print("The Python version doesn't the (built in) tkinter package. Please make sure you are using right Python interpreter. ( try: python3 -m packman gui )")
 
 
+#read and get_version is to print the PACKMAN version on the GUI
+def read(rel_path):
+    here = os.path.abspath(os.path.dirname(__file__))
+    with codecs.open(os.path.join(here, rel_path), 'r') as fp:
+        return fp.read()
+
+def get_version(rel_path):
+    '''To find out the version of the 
+    '''
+    for line in read(rel_path).splitlines():
+        if line.startswith('__version__'):
+            delim = '"' if '"' in line else "'"
+            return line.split(delim)[1]
+    else:
+        raise RuntimeError("Unable to find version string.")
+
 class Skeleton(tk.Tk):
     def __init__(self,frame_name):
         tk.Tk.__init__(self)
@@ -54,14 +73,14 @@ def __init__(self,frame_name):
 
         #Icon for windows
         try:
-            self.iconbitmap(False,'logo.ico')
+            self.iconbitmap(False, os.path.abspath(os.path.dirname(__file__)+'/logo.ico') )
         except:
             None
 
         self.show_frame(frame_name)
 
         try:
-            self.title_content = "PACKMAN GUI"
+            self.title_content = "PACKMAN "+get_version( os.path.abspath(os.path.dirname(__file__)).replace('bin','__init__.py') )+" GUI"
             self.title( self.title_content )
         except:
             self.title_content = "PACKMAN GUI"
@@ -533,6 +552,7 @@ def run_hdANM(self):
         tk.Label(pop_up1, text="Mode Number").grid(row=8,column=0,sticky=N+S+E+W, padx=10, pady=10 )
         tk.Label(pop_up1, text="Number of Frames").grid(row=9,column=0,sticky=N+S+E+W, padx=10, pady=10 )
         tk.Label(pop_up1, text="Movie Scale").grid(row=10,column=0,sticky=N+S+E+W, padx=10, pady=10 )
+        tk.Label(pop_up1, text="Project C-Alpha motions to all atoms").grid(row=11,column=0,sticky=N+S+E+W, padx=10, pady=10 )
 
         user_projection_choise = StringVar(pop_up1)
         user_projection_choise.set("Curvilinear")
@@ -545,13 +565,16 @@ def run_hdANM(self):
         n_frames.insert(END,10)
         n_scale  = tk.Entry(pop_up1)
         n_scale.grid(row=10,column=2,sticky=N+S+E+W, padx=10, pady=10 )
-        n_scale.insert(END,2)
-
-
-        self.save_movie_button = tk.Button(pop_up1,text='Save Movie',command=lambda : save_movie( n_mode, n_frames, n_scale, user_projection_choise.get() ) )
-        self.save_movie_button.grid(row=11,column=2,sticky=N+S+E+W, padx=10, pady=10 )
+        n_scale.insert(END,10)
+
+        #C-Alpha to all atom projection
+        ca_to_aa = BooleanVar()
+        n_ca_to_aa = tk.Checkbutton(pop_up1, variable=ca_to_aa, onvalue=True, offvalue= False)
+        n_ca_to_aa.grid(row=11,column=2,sticky=N+S+E+W, padx=10, pady=10 )
+
+        self.save_movie_button = tk.Button(pop_up1,text='Save Movie',command=lambda : save_movie( n_mode, n_frames, n_scale, user_projection_choise.get(), ca_to_aa ) )
+        self.save_movie_button.grid(row=12,column=2,sticky=N+S+E+W, padx=10, pady=10 )
 
-
         def save_matrix(choise):
             if(choise  == "hdANM Eigen Vectors"):
                 numpy.savetxt("hdANM Eigen Vectors.csv",Model.get_eigenvectors(),delimiter=",")
@@ -561,7 +584,7 @@ def save_matrix(choise):
                 showinfo('Notification','hdANM Eigen Values.csv file saved!')
             return True
 
-        def save_movie(n_mode, n_frames, n_scale, projection):
+        def save_movie(n_mode, n_frames, n_scale, projection, ca_to_aa):
             try:
                 n_mode = int(n_mode.get())
                 n_frames = int(n_frames.get())
@@ -570,9 +593,9 @@ def save_movie(n_mode, n_frames, n_scale, projection):
                 showinfo('Notification','Please check the parameters.\n\nMode Number (int)\nNumber of Frames (int)\nMovie Scale (float)')
 
             if(projection == "Curvilinear"):
-                Model.calculate_movie(n_mode, n=n_frames,scale=n_scale,extrapolation="curvilinear")
+                Model.calculate_movie(n_mode, n=n_frames,scale=n_scale,extrapolation="curvilinear",ca_to_aa=ca_to_aa.get())
             elif(projection =="Linear"):
-                Model.calculate_movie(n_mode, n=n_frames,scale=n_scale,extrapolation="linear")
+                Model.calculate_movie(n_mode, n=n_frames,scale=n_scale,extrapolation="linear",ca_to_aa=ca_to_aa.get())
             showinfo('Notification',str(n_mode)+'.pdb /.cif movie generated & saved.')
             return True
 

packman/bin/PACKMAN.py

@@ -61,6 +61,7 @@ def IO():
     hd_anm_io.add_argument("--scale", type=int, default=2, help='movie scale')
     hd_anm_io.add_argument("--frames", type=int, default=10, help='number of frames')
     hd_anm_io.add_argument("--modes", type=int, default=10, help='how many modes')
+    hd_anm_io.add_argument("--ca_to_aa", action=argparse.BooleanOptionalAction, type=bool, default=False, help='Project CA motion on all atoms.')
     hd_anm_io.add_argument("--make_tar", action='store_true', help='package output files into a tar.gz file')
 
     #Entropy