PDBAverage Python script added in the Preprocessing folder. The tool …

…has been also added to the GUI of the repository.

PELEAnalysis-Processing/PELEAnalysis_GUI

@@ -90,8 +90,9 @@ class App:
 		HStop.title("PDB Processor"); HStop.geometry("400x400")
 		self.bl = Label(HStop, text='', font = "helvetica 10"); self.bl.pack()
 		b1 = Button(HStop, text = "Preprocess", command = self.ppp, font = "helvetica 12")
-		b2 = Button(HStop, text = "Quit", command = HStop.destroy, font = "helvetica 12")
-		b1.pack(); b2.pack()
+		b2 = Button(HStop, text = "PDBAverage", command =self.pAv, font = "helvetica 12")
+		b3 = Button(HStop, text = "Quit", command = HStop.destroy, font = "helvetica 12")
+		b1.pack(); b2.pack(); b3.pack()
 
 	def ppp(self):
 		Input_filename = askopenfilename(title="The input PDB file", filetypes=[("PDB file",".pdb")])
@@ -105,6 +106,25 @@ class App:
 			output = subprocess.check_output("python Preprocessing/PDBProcessor4PELE.py -i {} -RN {}".format(Input_filename,Residue_name_list), shell = True)
 		self.bl['text'] = output.strip()
 
+	def pAv(self):
+		Input_filename = askopenfilename(title="The input PDB files (or trajectories)",
+			filetypes = [("PDB file",".pdb")])
+		Output_filename = askstring("The output filename", "Write down the output filename (Optional)")
+		BF = messagebox.askyesno("Do you want to calculate the B-factor for the average PDB file?")
+		if BF and Output_filename is not None:
+			output = subprocess.check_output("python Preprocessing/PDBAverage.py -i {} -o {} -B".format(
+				Input_filename, Output_filename), shell = True)
+		elif BF: 
+			output = subprocess.check_output("python Preprocessing/PDBAverage.py -i {} -B".format(
+				Input_filename), shell = True)
+		elif Output_filename is not None:
+			output = subprocess.check_output("python Preprocessing/PDBAverage.py -i {} -o {}".format(
+				Input_filename, Output_filename), shell = True)
+		else:
+			output = subprocess.check_output("python Preprocessing/PDBAverage.py -i {}".format(
+				Input_filename), shell = True)
+		self.bl['text'] = output.strip()
+
 	def PELE_scripting(self):
 		"""
 		Method that contains all the buttons and their actions

PELEAnalysis-Processing/Preprocessing/PDBAverage.py

@@ -0,0 +1,141 @@
+# -*- coding: utf-8 -*-
+
+# Global imports
+import os, time
+import glob
+import argparse as ap
+import numpy as np
+
+# Local imports
+from StorePDBFilenames import *
+
+# Script information
+__author__ = "Sergi Rodà Llordés"
+__version__ ="1.0"
+__maintainer__="Sergi Rodà Llordés"
+__email__="sergi.rodallordes@bsc.es"
+
+def parseArgs():
+    """
+    Parse arguments from command-line.
+
+    RETURNS
+    -------
+    PDBs: list
+            list of PDB files
+    """
+
+    parser = ap.ArgumentParser(description='Script used to average the coordinates from a set\
+        of PDB files')
+    optional = parser._action_groups.pop()
+    required = parser.add_argument_group('required arguments')
+    required.add_argument("-i", "--input", required=True, metavar="PDB FILE",
+                          type=str, nargs='*', help="path to PDB files")
+    optional.add_argument("-o", "--output", metavar="PDB FILE", type=str,
+                          help="the output PDB filename", default="average")
+    optional.add_argument("-B", "--B_factor", help="Output the calculated\
+        B_factor out of the average of coordinates", action="store_true")
+    parser._action_groups.append(optional)
+
+    args = parser.parse_args()
+
+    PDBs=storePDBfilenames(args.input, parser)
+
+    return PDBs, args.output, args.B_factor
+
+def getBFactor(atid_coords):
+    """
+    Estimate the B-factor from the root mean square fluctuation using
+    the average from the trajectories or list 
+    of PDB structurally aligned files.
+
+    RETURNS
+    -------
+    B_factor: string
+            Number of the temperature factor in string format
+    """
+
+    MSF = 0
+
+    for list_of_coordinates in atid_coords.values():
+        for coordinate in list_of_coordinates:
+            MSF += (abs(coordinate - np.average(list_of_coordinates)) ** 2) / len(list_of_coordinates)
+
+    B_factor = str(np.round(MSF * (8 * np.pi ** 2) / 3, 2))
+    if len(B_factor.split(".")[1])<2:
+        B_factor += (2 - len(B_factor.split(".")[1])) * "0"
+
+    return B_factor
+
+def averageStructure(PDBs, output, B_factor):
+    """
+    Function that calculates the average coordinates and 
+    outputs the "average" PDB file. The B-factor can be 
+    added to the atoms or not.
+
+    OUTPUT
+    ------
+    The average PDB file with the calculated B-factor or not
+    """
+
+    start = time.time()
+    D_av = {}
+    for PDB_filename in PDBs:
+        for atom in open(PDB_filename):
+            atid = atom[6:11]
+            if atid not in D_av:
+                D_av[atid] = {"x": [], "y": [], "z": []}
+            if atom[0:4] == "ATOM":
+                x = float(atom[30:38].strip(" "))
+                y = float(atom[38:46].strip(" "))
+                z = float(atom[46:54].strip(" "))
+                D_av[atid]["x"].append(x)
+                D_av[atid]["y"].append(y)
+                D_av[atid]["z"].append(z)
+
+
+    average_pdb = open(output + ".pdb","w")
+    for atom in open(PDBs[0]):
+        atid = atom[6:11]
+        if atom[0:4] == "ATOM":
+            x = str(round(np.average(D_av[atid]["x"]),3))
+            y = str(round(np.average(D_av[atid]["y"]),3))
+            z = str(round(np.average(D_av[atid]["z"]),3))
+            if len(x.split(".")[1])<3:
+                x += (3 - len(x.split(".")[1])) * "0"
+            if len(y.split(".")[1])<3:
+                y += (3 - len(y.split(".")[1])) * "0"
+            if len(z.split(".")[1])<3:
+                z += (3 - len(z.split(".")[1])) * "0"
+            if B_factor:
+                B_factor_value = getBFactor(D_av[atid])
+                average_pdb.write(atom[0:30]+"{:>8}".format(x)+"{:>8}".format(y)+"{:>8}".format(z)
+                +atom[54:60]+"{:>6}".format(B_factor_value)+atom[66:])
+            else:
+                average_pdb.write(atom[0:30]+"{:>8}".format(x)+"{:>8}".format(y)+"{:>8}".format(z)
+                +atom[54:])
+        elif atom[0:3] == "TER" or atom[0:3] == "END" or atom[0:6] == "HETATM":
+            average_pdb.write("TER")
+            break
+        else:
+            continue
+
+    end = time.time()
+    print("\nThe main code needed {} seconds to compute the average PDB file \n".format(end-start))
+
+def main():
+    """
+    Main function
+
+    It is called when this script is the main program called by the interpreter.
+    """
+
+    # Parse command-line arguments
+    PDBs, output, B_factor = parseArgs()
+
+    # Average PDB structure from input PDB files
+    averageStructure(PDBs, output, B_factor)
+
+if __name__=="__main__":
+	"""Call the main function"""
+	main()

PELEAnalysis-Processing/Preprocessing/PDBProcessor4PELE.py

@@ -5,6 +5,9 @@
 import glob
 import argparse as ap
 
+# Local imports
+from StorePDBFilenames import *
+
 # Script information
 __author__ = "Sergi Rodà"
 __license__ = "MIT"
@@ -20,37 +23,6 @@
  process them to be able to be used
  in a PELE simulation"""
 
-def storePDBfilenames(PDBs_to_parse, parser):
-    """
-    It identifies the PDB files to add to the PDBProcessor4PELE tool
-
-    PARAMETERS
-    ----------
-    PDBs_to_parse : list of strings
-                       all the PDB files that want to be added to the analysis
-
-    RETURNS
-    -------
-    parsed_data : list of PDB filenames (strings)
-    """
-
-    PDBs = []
-
-    for PDB_list in PDBs_to_parse:
-        PDB_found = glob.glob(PDB_list)
-        if len(PDB_found) == 0:
-            print("Warning: path to PDB file \'" +
-                  "{}".format(PDB_list) + "\' not found.")
-        for PDB in PDB_found:
-            PDBs.append(PDB)
-
-    if len(PDBs) == 0:
-        print("Error: list of PDB files is empty.")
-        parser.print_help()
-        exit(1)
-
-    return PDBs
-
 def parseArgs():
     """
     Parse arguments from command-line

PELEAnalysis-Processing/Preprocessing/StorePDBFilenames.py

@@ -0,0 +1,44 @@
+# -*- coding: utf-8 -*-
+
+# Global imports
+from __future__ import unicode_literals
+import os
+import glob
+
+# Script information
+__author__ = "Sergi Rodà"
+__license__ = "MIT"
+__version__ = "1.0.1"
+__maintainer__ = "Sergi Rodà"
+__email__ = "sergi.rodallordes@bsc.es"
+
+def storePDBfilenames(PDBs_to_parse, parser):
+    """
+    It identifies the PDB files to add to the PDBProcessor4PELE tool
+
+    PARAMETERS
+    ----------
+    PDBs_to_parse : list of strings
+                       all the PDB files that want to be added to the analysis
+
+    RETURNS
+    -------
+    parsed_data : list of PDB filenames (strings)
+    """
+
+    PDBs = []
+
+    for PDB_list in PDBs_to_parse:
+        PDB_found = glob.glob(PDB_list)
+        if len(PDB_found) == 0:
+            print("Warning: path to PDB file \'" +
+                  "{}".format(PDB_list) + "\' not found.")
+        for PDB in PDB_found:
+            PDBs.append(PDB)
+
+    if len(PDBs) == 0:
+        print("Error: list of PDB files is empty.")
+        parser.print_help()
+        exit(1)
+
+    return PDBs