From 2d62238bc15d43eef005e1a9491c632b587ae79b Mon Sep 17 00:00:00 2001 From: jgilaber Date: Tue, 27 Oct 2020 11:19:57 +0100 Subject: [PATCH 01/27] Add helper methods to Dihedral object --- offpele/topology/topology.py | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/offpele/topology/topology.py b/offpele/topology/topology.py index 47ec7c04..81921161 100644 --- a/offpele/topology/topology.py +++ b/offpele/topology/topology.py @@ -465,6 +465,19 @@ def set_atom4_idx(self, index): """ self._atom4_idx = index + def get_bonds(self): + """ + It returns a list of tuples containing the atom indices of the bonds + that make the dihedral + + Returns + ------- + bonds : list + List of the bonds involved in this Dihedral object, identified by + the indexes of the atoms + """ + return [(self._atom1_idx, self._atom2_idx), (self.atom2_idx, self._atom3_idx), (self.atom3_idx, self.atom4_idx)] + def plot(self): """ It plots this Dihedral as a function of phi angle. @@ -540,6 +553,18 @@ def atom4_idx(self): """ return self._atom4_idx + @property + def atoms(self): + """ + Dihedral's atom indexes + + Returns + ------- + atoms_idx: tuple + The indexes of all atoms involved in this Dihedral object + """ + return (self._atom1_idx, self._atom2_idx, self._atom3_idx, self._atom4_idx) + @property def periodicity(self): """ From fff396985d58c48fc886ebe9ea54e320db45fe70 Mon Sep 17 00:00:00 2001 From: jgilaber Date: Tue, 27 Oct 2020 17:28:14 +0100 Subject: [PATCH 02/27] Add options to main to run dihedral calculations --- offpele/main.py | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/offpele/main.py b/offpele/main.py index b3d8b115..fe742479 100644 --- a/offpele/main.py +++ b/offpele/main.py @@ -54,6 +54,9 @@ def parse_args(args): parser.add_argument("-o", "--output", metavar="PATH", help="Output path. Default is the current working " + "directory") + parser.add_argument("--dihedrals_info_path", default=None, type=str, + help="Path to the folder containing the BCE output" + + " used to collect dihedral angles for PELE") parser.add_argument('--with_solvent', dest='with_solvent', help="Generate solvent parameters for OBC", action='store_true') @@ -93,7 +96,8 @@ def run_offpele(pdb_file, forcefield=DEFAULT_OFF_FORCEFIELD, resolution=DEFAULT_RESOLUTION, charge_method=DEFAULT_CHARGE_METHOD, exclude_terminal_rotamers=True, - output=None, with_solvent=False, as_datalocal=False): + output=None, with_solvent=False, as_datalocal=False, + dihedral_path=None): """ It runs offpele. @@ -139,6 +143,7 @@ def run_offpele(pdb_file, forcefield=DEFAULT_OFF_FORCEFIELD, from offpele.topology import Molecule from offpele.template import Impact from offpele.solvent import OBC2 + from offpele.BCEDihedrals import BCEDihedrals if not output: output = os.getcwd() @@ -160,6 +165,11 @@ def run_offpele(pdb_file, forcefield=DEFAULT_OFF_FORCEFIELD, solvent = OBC2(molecule) solvent.to_json_file(output_handler.get_solvent_template_path()) + if dihedral_path is not None: + dihedrals = BCEDihedrals(molecule, dihedral_path) + dihedrals.calculate() + dihedrals.save(output_handler.get_dihedral_library_path()) + log.info(' - All files were generated successfully') log.info('-' * 60) @@ -200,7 +210,8 @@ def main(args): run_offpele(args.pdb_file, args.forcefield, args.resolution, args.charge_method, exclude_terminal_rotamers, - args.output, args.with_solvent, args.as_datalocal) + args.output, args.with_solvent, args.as_datalocal, + dihedral_path=args.dihedrals_info_path) if __name__ == '__main__': From f0e2a72b078e4a08ddbd837acabad57ac906e76f Mon Sep 17 00:00:00 2001 From: jgilaber Date: Tue, 27 Oct 2020 17:29:05 +0100 Subject: [PATCH 03/27] Add default path to dihedral library file --- offpele/utils/output.py | 32 +++++++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/offpele/utils/output.py b/offpele/utils/output.py index f264c6d7..b36ababe 100644 --- a/offpele/utils/output.py +++ b/offpele/utils/output.py @@ -18,7 +18,8 @@ class OutputPathHandler(object): OPLS_IMPACT_TEMPLATE_PATH = 'DataLocal/Templates/OPLS2005/HeteroAtoms/' ROTAMER_LIBRARY_PATH = 'DataLocal/LigandRotamerLibs/' SOLVENT_TEMPLATE_PATH = 'DataLocal/OBC/' - FILE_TYPES = ['impact template', 'rotamer library', 'solvent template'] + DIHEDRALS_LIBRARY_PATH = 'DataLocal/Dihedrals/' + FILE_TYPES = ['impact template', 'rotamer library', 'solvent template', 'dihedral library'] def __init__(self, molecule, as_datalocal=False, output_path=None): """ @@ -84,6 +85,9 @@ def get_path(self, file_type, create_missing_folders=True): if file_type.lower() == 'solvent template': return self.get_solvent_template_path(create_missing_folders) + if file_type.lower() == 'dihedral library': + return self.get_dihedral_library_path(create_missing_folders) + def get_impact_template_path(self, create_missing_folders=True): """ It returns the path for an Impact template file. @@ -176,6 +180,32 @@ def get_solvent_template_path(self, create_missing_folders=True): return os.path.join(path, file_name) + def get_dihedral_library_path(self, create_missing_folders=True): + """ + It returns the path for a dihedral library file. + + Parameters + ---------- + create_missing_folders : bool + Whether to create missing folders or not. Default is True + + Returns + ------- + file_path : str + The path for a dihedral library file + """ + file_name = self._molecule.tag.upper() + '.dihedral' + + if self.as_datalocal: + path = os.path.join(self.output_path, self.DIHEDRALS_LIBRARY_PATH) + else: + path = self.output_path + + if create_missing_folders: + create_path(path) + + return os.path.join(path, file_name) + @property def as_datalocal(self): """ From 6717f0cf91b6ec247f6a62777f0d6eebc2159c08 Mon Sep 17 00:00:00 2001 From: jgilaber Date: Tue, 27 Oct 2020 17:29:28 +0100 Subject: [PATCH 04/27] Add module to calculate dihedral angles --- offpele/BCEDihedrals/BCEDihedrals.py | 118 +++++++++++++++++++++++++++ offpele/BCEDihedrals/__init__.py | 1 + 2 files changed, 119 insertions(+) create mode 100644 offpele/BCEDihedrals/BCEDihedrals.py create mode 100644 offpele/BCEDihedrals/__init__.py diff --git a/offpele/BCEDihedrals/BCEDihedrals.py b/offpele/BCEDihedrals/BCEDihedrals.py new file mode 100644 index 00000000..3cc392a0 --- /dev/null +++ b/offpele/BCEDihedrals/BCEDihedrals.py @@ -0,0 +1,118 @@ +""" +This module contains classes and functions involved in the +creation of a library of dihedral angles from the output of the +BCE server +""" +import os +import glob +from rdkit.Chem import rdMolTransforms +from offpele.utils import Logger +from offpele.topology import molecule + +class BCEDihedrals(object): + """ + A class to produce a library of dihedral angles from the output + of the BCE server + """ + def __init__(self, molecule_obj, bce_output_path): + """ + Initializes a BCEDihedrals object + + Parameters + ---------- + molecule_obj : An offpele.topology.Molecule + A Molecule object that contains the ligand's information + + bce_output_path: str + Path where the output from the BCE server is stored + """ + self._molecule = molecule_obj + self.bce_path = bce_output_path + self.dihedral_library = {} + + def calculate(self): + """ + Calculate dihedrals library from the bce output + """ + logger = Logger() + logger.info(' - Calculating dihedral library') + # here we rely on the output of the bce server to not change its + # structrue in the future + dihedral_file = glob.glob(os.path.join(self.bce_path, "DIHEDRALS", "*_new.ndx"))[0] + dihedrals = read_dihedral_info_file(dihedral_file) + clusters = glob.glob(os.path.join(self.bce_path, "CLUSTERS", "CL*", "cluster*.min.imaged.pdb")) + for cluster in clusters: + self.calculate_cluster_angles(cluster, dihedrals) + + + def calculate_cluster_angles(self, cluster_pdb, dihedral_list): + """ + Calculate dihedral angles from pdb + + Parameters + ---------- + cluster_pdb: str + Path to the cluster representative conformation + + dihedral_list: list + List of the tuples containing the atoms that form the dihedrals + """ + pdb_dihedrals = [] + # use the input molecule as template since the cluster structures + # probably will not have proper stereochemistry + mol = molecule.Molecule(cluster_pdb, connectivity_template=self._molecule.rdkit_molecule) + # we use substructure matching to ensure that the indices in the + # clusters pdb and the input ligand are the same + rename_dict = {i: x for i, x in enumerate(self._molecule.rdkit_molecule.GetSubstructMatch(mol.rdkit_molecule))} + for dihedral in dihedral_list: + names = [self._molecule.atoms[rename_dict[atom]].PDB_name for atom in dihedral] + angle = rdMolTransforms.GetDihedralRad(mol.rdkit_molecule.GetConformer(), *dihedral) + pdb_dihedrals.append(names+[angle]) + self.dihedral_library[cluster_pdb] = pdb_dihedrals + + def save(self, output_path): + """ + Save the dihedral library + + Parameters + ---------- + output_path: str + Where to save the dihedral library + """ + with open(output_path, "w") as fw: + for dihedral_file, dihedral_values in self.dihedral_library.items(): + fw.write("BEGINDIHEDRALS\n") + fw.write("File: {:s}\n".format(dihedral_file)) + fw.write("Atom 1 Atom 2 Atom3 Atom 4 Angle \n") + for dihedral in dihedral_values: + fw.write("{:s} {:s} {:s} {:s} {:5f}\n".format(*dihedral)) + fw.write("ENDDIHEDRALS\n") + fw.write("END") + +def read_dihedral_info_file(file_name): + """ + Read an ndx file containing the classification of all dihedrals and return + the flexible ones + + Parameters + ---------- + file_name: str + Path of the dihedrals file + """ + read_dihedral = False + dihedrals = [] + with open(file_name) as f: + for line in f: + # reached the block of all dihedrals, we can stop reading + if line.startswith("[ ALL ]"): + return dihedrals + if line.startswith("["): + name = line.strip("\n[] ") + read_dihedral = name.startswith("F") + else: + if read_dihedral: + # the indices of the bce server are 1-based, while our + # indices are 0-based + dihedrals.append([int(x)-1 for x in line.rstrip().split()]) + read_dihedral = False + return dihedrals diff --git a/offpele/BCEDihedrals/__init__.py b/offpele/BCEDihedrals/__init__.py new file mode 100644 index 00000000..1d38f1af --- /dev/null +++ b/offpele/BCEDihedrals/__init__.py @@ -0,0 +1 @@ +from .BCEDihedrals import BCEDihedrals From 93fee125044f03804d420c531de741f97094ff8d Mon Sep 17 00:00:00 2001 From: jgilaber Date: Wed, 28 Oct 2020 17:48:34 +0100 Subject: [PATCH 05/27] Move rdkit dependencies into its toolkit --- offpele/BCEDihedrals/BCEDihedrals.py | 15 ++++++---- offpele/utils/toolkits.py | 41 ++++++++++++++++++++++++++++ 2 files changed, 50 insertions(+), 6 deletions(-) diff --git a/offpele/BCEDihedrals/BCEDihedrals.py b/offpele/BCEDihedrals/BCEDihedrals.py index 3cc392a0..64023bee 100644 --- a/offpele/BCEDihedrals/BCEDihedrals.py +++ b/offpele/BCEDihedrals/BCEDihedrals.py @@ -5,9 +5,10 @@ """ import os import glob -from rdkit.Chem import rdMolTransforms +import offpele from offpele.utils import Logger from offpele.topology import molecule +from offpele.utils.toolkits import RDKitToolkitWrapper class BCEDihedrals(object): """ @@ -57,16 +58,17 @@ def calculate_cluster_angles(self, cluster_pdb, dihedral_list): dihedral_list: list List of the tuples containing the atoms that form the dihedrals """ + rdkit_wrapper = RDKitToolkitWrapper() pdb_dihedrals = [] # use the input molecule as template since the cluster structures # probably will not have proper stereochemistry mol = molecule.Molecule(cluster_pdb, connectivity_template=self._molecule.rdkit_molecule) # we use substructure matching to ensure that the indices in the # clusters pdb and the input ligand are the same - rename_dict = {i: x for i, x in enumerate(self._molecule.rdkit_molecule.GetSubstructMatch(mol.rdkit_molecule))} + rename_dict = {i: x for i, x in enumerate(rdkit_wrapper.get_substruct_match(self._molecule, mol))} for dihedral in dihedral_list: names = [self._molecule.atoms[rename_dict[atom]].PDB_name for atom in dihedral] - angle = rdMolTransforms.GetDihedralRad(mol.rdkit_molecule.GetConformer(), *dihedral) + angle = rdkit_wrapper.get_dihedral(mol, *dihedral, units="degrees") pdb_dihedrals.append(names+[angle]) self.dihedral_library[cluster_pdb] = pdb_dihedrals @@ -80,10 +82,11 @@ def save(self, output_path): Where to save the dihedral library """ with open(output_path, "w") as fw: + fw.write("* DIHEDRAL LIBRARY FILE\n") + fw.write('* File generated with offpele-{}\n'.format( offpele.__version__)) for dihedral_file, dihedral_values in self.dihedral_library.items(): - fw.write("BEGINDIHEDRALS\n") - fw.write("File: {:s}\n".format(dihedral_file)) - fw.write("Atom 1 Atom 2 Atom3 Atom 4 Angle \n") + fw.write("* File: {:s}\n".format(dihedral_file)) + fw.write("{:s} {:d} {:d}\n".format(self._molecule.tag.upper(), len(dihedral_values), len(self.dihedral_library))) for dihedral in dihedral_values: fw.write("{:s} {:s} {:s} {:s} {:5f}\n".format(*dihedral)) fw.write("ENDDIHEDRALS\n") diff --git a/offpele/utils/toolkits.py b/offpele/utils/toolkits.py index f1015565..51d77ea0 100644 --- a/offpele/utils/toolkits.py +++ b/offpele/utils/toolkits.py @@ -264,6 +264,47 @@ def get_atom_names(self, molecule): return atom_names + def get_dihedral(self, mol, atom1, atom2, atom3, atom4, units="radians"): + """ + It calculates the value of the dihedral angle in the specified units + (default radians) + + Parameters + ---------- + molecule : an offpele.topology.Molecule + The offpele's Molecule object + atom1 : int + Index of the first atom in the dihedral + atom2 : int + Index of the second atom in the dihedral + atom3 : int + Index of the third atom in the dihedral + atom4 : int + Index of the fourth atom in the dihedral + units : str + The units in which to calculate the angle (default is radians, can + be radians or degrees) + """ + from rdkit.Chem import rdMolTransforms + if units == "degrees": + angle = rdMolTransforms.GetDihedralDeg(mol.rdkit_molecule.GetConformer(), atom1, atom2, atom3, atom4) + else: + angle = rdMolTransforms.GetDihedralRad(mol.rdkit_molecule.GetConformer(), atom1, atom2, atom3, atom4) + return angle + + def get_substruct_match(self, mol1, mol2): + """ + It returns the atoms in mol2 that match those of mol1 + + Parameters + ---------- + mol1 : an offpele.topology.Molecule + Molecule to match atoms from + mol2 : an offpele.topology.Molecule + Molecule with the atoms to match + """ + return mol1.rdkit_molecule.GetSubstructMatch(mol2.rdkit_molecule) + def to_pdb_file(self, molecule, path): """ It writes the RDKit molecule to a PDB file. From d7ebbe1e05b95a651e01d57a8390aadc95ef4b3d Mon Sep 17 00:00:00 2001 From: jgilaber Date: Fri, 30 Oct 2020 17:35:15 +0100 Subject: [PATCH 06/27] Change units of dihedral calculation to rad --- offpele/BCEDihedrals/BCEDihedrals.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/offpele/BCEDihedrals/BCEDihedrals.py b/offpele/BCEDihedrals/BCEDihedrals.py index 64023bee..a0a29fa7 100644 --- a/offpele/BCEDihedrals/BCEDihedrals.py +++ b/offpele/BCEDihedrals/BCEDihedrals.py @@ -68,7 +68,7 @@ def calculate_cluster_angles(self, cluster_pdb, dihedral_list): rename_dict = {i: x for i, x in enumerate(rdkit_wrapper.get_substruct_match(self._molecule, mol))} for dihedral in dihedral_list: names = [self._molecule.atoms[rename_dict[atom]].PDB_name for atom in dihedral] - angle = rdkit_wrapper.get_dihedral(mol, *dihedral, units="degrees") + angle = rdkit_wrapper.get_dihedral(mol, *dihedral, units="radians") pdb_dihedrals.append(names+[angle]) self.dihedral_library[cluster_pdb] = pdb_dihedrals From 4720fe86fa63b2b76fe149f9048c7b7478f5949d Mon Sep 17 00:00:00 2001 From: cescgina Date: Fri, 27 Nov 2020 20:12:36 +0100 Subject: [PATCH 07/27] Remove duplicated method from merge --- peleffy/topology/topology.py | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/peleffy/topology/topology.py b/peleffy/topology/topology.py index 7c0e2fc1..c94729d4 100644 --- a/peleffy/topology/topology.py +++ b/peleffy/topology/topology.py @@ -374,17 +374,6 @@ def atoms(self): return self._atoms @property - def atoms(self): - """ - Dihedral's atom indexes - - Returns - ------- - atoms_idx: tuple - The indexes of all atoms involved in this Dihedral object - """ - return (self._atom1_idx, self._atom2_idx, self._atom3_idx, self._atom4_idx) - def bonds(self): """ The list of bonds in the topology. From 8a3b62a82515423b4100800882740f0c2d57c81a Mon Sep 17 00:00:00 2001 From: cescgina Date: Fri, 27 Nov 2020 20:14:10 +0100 Subject: [PATCH 08/27] Fix discrepancy in name of OPLS charges option --- peleffy/main.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/peleffy/main.py b/peleffy/main.py index d116c32a..af607378 100644 --- a/peleffy/main.py +++ b/peleffy/main.py @@ -19,7 +19,7 @@ DEFAULT_OFF_FORCEFIELD = 'openff_unconstrained-1.2.1.offxml' DEFAULT_RESOLUTION = int(30) DEFAULT_CHARGE_METHOD = 'am1bcc' -AVAILABLE_CHARGE_METHODS = ['am1bcc', 'gasteiger', 'OPLS'] +AVAILABLE_CHARGE_METHODS = ['am1bcc', 'gasteiger', 'OPLS2005'] IMPACT_TEMPLATE_PATH = 'DataLocal/Templates/OFF/Parsley/HeteroAtoms/' ROTAMER_LIBRARY_PATH = 'DataLocal/LigandRotamerLibs/' SOLVENT_TEMPLATE_PATH = 'DataLocal/OBC/' From 23bdca88786eb3d052e854084998524f431a4ff0 Mon Sep 17 00:00:00 2001 From: cescgina Date: Fri, 27 Nov 2020 20:30:19 +0100 Subject: [PATCH 09/27] Udpate format of Dihedral library due to new changes --- peleffy/BCEDihedrals/BCEDihedrals.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/peleffy/BCEDihedrals/BCEDihedrals.py b/peleffy/BCEDihedrals/BCEDihedrals.py index c4aeda37..e6ea6372 100644 --- a/peleffy/BCEDihedrals/BCEDihedrals.py +++ b/peleffy/BCEDihedrals/BCEDihedrals.py @@ -67,7 +67,7 @@ def calculate_cluster_angles(self, cluster_pdb, dihedral_list): # clusters pdb and the input ligand are the same rename_dict = {i: x for i, x in enumerate(rdkit_wrapper.get_substruct_match(self._molecule, mol))} for dihedral in dihedral_list: - names = [self._molecule.atoms[rename_dict[atom]].PDB_name for atom in dihedral] + names = [self._molecule.get_pdb_atom_names()[rename_dict[atom]] for atom in dihedral] angle = rdkit_wrapper.get_dihedral(mol, *dihedral, units="radians") pdb_dihedrals.append(names+[angle]) self.dihedral_library[cluster_pdb] = pdb_dihedrals @@ -88,7 +88,7 @@ def save(self, output_path): fw.write("* File: {:s}\n".format(dihedral_file)) fw.write("{:s} {:d} {:d}\n".format(self._molecule.tag.upper(), len(dihedral_values), len(self.dihedral_library))) for dihedral in dihedral_values: - fw.write("{:s} {:s} {:s} {:s} {:5f}\n".format(*dihedral)) + fw.write("{:s}/{:s}/{:s}/{:s}/{:5f}\n".format(*dihedral)) fw.write("ENDDIHEDRALS\n") fw.write("END") From 588e8ec50250a869ab0a8b281a43a5f81fdd97d6 Mon Sep 17 00:00:00 2001 From: cescgina Date: Fri, 27 Nov 2020 20:50:19 +0100 Subject: [PATCH 10/27] Revert changes in dihedral library format By using the topology instead of the molecule we can retain the same format as before --- peleffy/BCEDihedrals/BCEDihedrals.py | 14 ++++++++------ peleffy/main.py | 2 +- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/peleffy/BCEDihedrals/BCEDihedrals.py b/peleffy/BCEDihedrals/BCEDihedrals.py index e6ea6372..a8cb9226 100644 --- a/peleffy/BCEDihedrals/BCEDihedrals.py +++ b/peleffy/BCEDihedrals/BCEDihedrals.py @@ -15,19 +15,20 @@ class BCEDihedrals(object): A class to produce a library of dihedral angles from the output of the BCE server """ - def __init__(self, molecule_obj, bce_output_path): + def __init__(self, topology_obj, bce_output_path): """ Initializes a BCEDihedrals object Parameters ---------- - molecule_obj : An peleffy.topology.Molecule - A Molecule object that contains the ligand's information + molecule_obj : An peleffy.topology.Topology + A Topology object that contains the ligand's information bce_output_path: str Path where the output from the BCE server is stored """ - self._molecule = molecule_obj + self._topology = topology_obj + self._molecule = topology_obj.molecule self.bce_path = bce_output_path self.dihedral_library = {} @@ -67,7 +68,8 @@ def calculate_cluster_angles(self, cluster_pdb, dihedral_list): # clusters pdb and the input ligand are the same rename_dict = {i: x for i, x in enumerate(rdkit_wrapper.get_substruct_match(self._molecule, mol))} for dihedral in dihedral_list: - names = [self._molecule.get_pdb_atom_names()[rename_dict[atom]] for atom in dihedral] + # names = [self._molecule.get_pdb_atom_names()[rename_dict[atom]] for atom in dihedral] + names = [self._topology.atoms[rename_dict[atom]].PDB_name for atom in dihedral] angle = rdkit_wrapper.get_dihedral(mol, *dihedral, units="radians") pdb_dihedrals.append(names+[angle]) self.dihedral_library[cluster_pdb] = pdb_dihedrals @@ -88,7 +90,7 @@ def save(self, output_path): fw.write("* File: {:s}\n".format(dihedral_file)) fw.write("{:s} {:d} {:d}\n".format(self._molecule.tag.upper(), len(dihedral_values), len(self.dihedral_library))) for dihedral in dihedral_values: - fw.write("{:s}/{:s}/{:s}/{:s}/{:5f}\n".format(*dihedral)) + fw.write("{:s} {:s} {:s} {:s} {:5f}\n".format(*dihedral)) fw.write("ENDDIHEDRALS\n") fw.write("END") diff --git a/peleffy/main.py b/peleffy/main.py index af607378..ce17b56d 100644 --- a/peleffy/main.py +++ b/peleffy/main.py @@ -183,7 +183,7 @@ def run_peleffy(pdb_file, solvent.to_file(output_handler.get_solvent_template_path()) if dihedral_path is not None: - dihedrals = BCEDihedrals(molecule, dihedral_path) + dihedrals = BCEDihedrals(topology, dihedral_path) dihedrals.calculate() dihedrals.save(output_handler.get_dihedral_library_path()) From c2c666a17025992daca999311c5355ba6d90532a Mon Sep 17 00:00:00 2001 From: cescgina Date: Fri, 27 Nov 2020 20:56:55 +0100 Subject: [PATCH 11/27] Sort the output of glob to get consistent output This will produce the exact same dihedral library in all machines --- peleffy/BCEDihedrals/BCEDihedrals.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/peleffy/BCEDihedrals/BCEDihedrals.py b/peleffy/BCEDihedrals/BCEDihedrals.py index a8cb9226..eac7e9e6 100644 --- a/peleffy/BCEDihedrals/BCEDihedrals.py +++ b/peleffy/BCEDihedrals/BCEDihedrals.py @@ -42,7 +42,7 @@ def calculate(self): # structrue in the future dihedral_file = glob.glob(os.path.join(self.bce_path, "DIHEDRALS", "*_new.ndx"))[0] dihedrals = read_dihedral_info_file(dihedral_file) - clusters = glob.glob(os.path.join(self.bce_path, "CLUSTERS", "CL*", "cluster*.min.imaged.pdb")) + clusters = sorted(glob.glob(os.path.join(self.bce_path, "CLUSTERS", "CL*", "cluster*.min.imaged.pdb"))) for cluster in clusters: self.calculate_cluster_angles(cluster, dihedrals) @@ -68,7 +68,6 @@ def calculate_cluster_angles(self, cluster_pdb, dihedral_list): # clusters pdb and the input ligand are the same rename_dict = {i: x for i, x in enumerate(rdkit_wrapper.get_substruct_match(self._molecule, mol))} for dihedral in dihedral_list: - # names = [self._molecule.get_pdb_atom_names()[rename_dict[atom]] for atom in dihedral] names = [self._topology.atoms[rename_dict[atom]].PDB_name for atom in dihedral] angle = rdkit_wrapper.get_dihedral(mol, *dihedral, units="radians") pdb_dihedrals.append(names+[angle]) From 3b6bea39ebff9e0afc67f7c72a4216f69e494cf3 Mon Sep 17 00:00:00 2001 From: cescgina Date: Sat, 28 Nov 2020 11:26:50 +0100 Subject: [PATCH 12/27] Expose charge calculator types to main for consitency --- peleffy/main.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/peleffy/main.py b/peleffy/main.py index 3f51d3b0..bab4fe11 100644 --- a/peleffy/main.py +++ b/peleffy/main.py @@ -14,12 +14,13 @@ import peleffy from peleffy.utils import Logger, OutputPathHandler +from peleffy.forcefield.selectors import ChargeCalculatorSelector DEFAULT_OFF_FORCEFIELD = 'openff_unconstrained-1.3.0.offxml' DEFAULT_RESOLUTION = int(30) DEFAULT_CHARGE_METHOD = 'am1bcc' -AVAILABLE_CHARGE_METHODS = ['am1bcc', 'gasteiger', 'OPLS2005'] +AVAILABLE_CHARGE_METHODS = ChargeCalculatorSelector()._AVAILABLE_TYPES.keys() IMPACT_TEMPLATE_PATH = 'DataLocal/Templates/OFF/Parsley/HeteroAtoms/' ROTAMER_LIBRARY_PATH = 'DataLocal/LigandRotamerLibs/' SOLVENT_TEMPLATE_PATH = 'DataLocal/OBC/' From 4e1aff9bb9f967c95cba1adacfd7bab50f08e806 Mon Sep 17 00:00:00 2001 From: jgilaber Date: Mon, 7 Dec 2020 15:12:22 +0100 Subject: [PATCH 13/27] Add function to calculate RMSD in RDKIT toolkit --- peleffy/utils/toolkits.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/peleffy/utils/toolkits.py b/peleffy/utils/toolkits.py index 0d1b3cc7..fd82efd9 100644 --- a/peleffy/utils/toolkits.py +++ b/peleffy/utils/toolkits.py @@ -568,6 +568,26 @@ def get_2D_representation(self, molecule): AllChem.Compute2DCoords(representation_2D) return representation_2D + def get_rmsd(self, molecule, molecule_2): + """ + It returns the RMSD between two RDKit molecules. + + Parameters + ---------- + molecule : an peleffy.topology.Molecule + The peleffy's Molecule object + molecule_2 : an peleffy.topology.Molecule + The peleffy's Molecule object + + Returns + ------- + rmsd_value : float + RMSD between two RDKit molecules + """ + from rdkit.Chem import rdMolAlign + return rdMolAlign.AlignMol(molecule.rdkit_molecule, + molecule_2.rdkit_molecule) + def draw_molecule(self, representation, atom_indexes=list(), radii_dict=dict(), atom_color_dict=dict(), bond_indexes=list(), bond_color_dict=dict()): From 24a03867b493451cccce263b2b1191b532b41ef1 Mon Sep 17 00:00:00 2001 From: jgilaber Date: Mon, 7 Dec 2020 15:23:52 +0100 Subject: [PATCH 14/27] Modify logger level inside BCEDihedrals --- peleffy/main.py | 4 +++- peleffy/utils/utils.py | 24 ++++++++++++++++++++++++ 2 files changed, 27 insertions(+), 1 deletion(-) diff --git a/peleffy/main.py b/peleffy/main.py index bab4fe11..7d7db030 100644 --- a/peleffy/main.py +++ b/peleffy/main.py @@ -191,9 +191,11 @@ def run_peleffy(pdb_file, solvent.to_file(output_handler.get_solvent_template_path()) if dihedral_path is not None: - dihedrals = BCEDihedrals(topology, dihedral_path) + previous_level = log.get_level() + dihedrals = BCEDihedrals(topology, dihedral_path, mode=dihedral_mode) dihedrals.calculate() dihedrals.save(output_handler.get_dihedral_library_path()) + log.set_level(previous_level) log.info(' - All files were generated successfully:') log.info(' - {}'.format(output_handler.get_rotamer_library_path())) diff --git a/peleffy/utils/utils.py b/peleffy/utils/utils.py index de68bcd2..bac41d2b 100644 --- a/peleffy/utils/utils.py +++ b/peleffy/utils/utils.py @@ -279,6 +279,30 @@ def set_level(self, level): for handler in self._logger.handlers: handler.setLevel(logging_level) + def get_level(self): + """ + It gets the logging level. + + Returns + ------- + level : str + The logging level to set. One of [DEBUG, INFO, WARNING, ERROR, + CRITICAL] + """ + import logging + level = self._logger.level + + if level == logging.DEBUG: + return 'DEBUG' + if level == logging.INFO: + return 'INFO' + if level == logging.WARNING: + return 'WARNING' + if level == logging.ERROR: + return 'ERROR' + if level == logging.CRITICAL: + return 'CRITICAL' + def debug(self, *messages): """ It pulls a debug message. From 9122c8596828b7e8bb2a0a61dc65620e2e845659 Mon Sep 17 00:00:00 2001 From: jgilaber Date: Mon, 7 Dec 2020 15:24:12 +0100 Subject: [PATCH 15/27] Add new mode in BCEDihedrals In this new mode, the different clusters are ordered so that the jump between two subsequent clusters involves the minimum structural difference. This is done heuristically, and it is not guaranteed to be optimal --- peleffy/BCEDihedrals/BCEDihedrals.py | 131 +++++++++++++++++++++++++-- peleffy/main.py | 10 +- 2 files changed, 132 insertions(+), 9 deletions(-) diff --git a/peleffy/BCEDihedrals/BCEDihedrals.py b/peleffy/BCEDihedrals/BCEDihedrals.py index eac7e9e6..e2223b13 100644 --- a/peleffy/BCEDihedrals/BCEDihedrals.py +++ b/peleffy/BCEDihedrals/BCEDihedrals.py @@ -5,6 +5,8 @@ """ import os import glob +import numpy as np +import networkx as nx import peleffy from peleffy.utils import Logger from peleffy.topology import molecule @@ -15,22 +17,26 @@ class BCEDihedrals(object): A class to produce a library of dihedral angles from the output of the BCE server """ - def __init__(self, topology_obj, bce_output_path): + def __init__(self, topology_obj, bce_output_path, mode="all_dihedrals", verbose=False): """ Initializes a BCEDihedrals object Parameters ---------- - molecule_obj : An peleffy.topology.Topology + topology_obj : An peleffy.topology.Topology A Topology object that contains the ligand's information - bce_output_path: str Path where the output from the BCE server is stored + mode: str + Whether to extract all dihedrals or only those marked as flexible + """ self._topology = topology_obj self._molecule = topology_obj.molecule self.bce_path = bce_output_path self.dihedral_library = {} + self.mode = mode + self.verbose = verbose def calculate(self): """ @@ -38,6 +44,23 @@ def calculate(self): """ logger = Logger() logger.info(' - Calculating dihedral library') + logger.info(' - with mode: {}'.format(self.mode)) + if not self.verbose: + logger.set_level('WARNING') + if self.mode == "all_dihedrals": + self._calculate_all_dihedrals() + elif self.mode == "flexible_dihedrals": + self._calculate_flexible_dihedrals() + + def _calculate_all_dihedrals(self): + clusters = sorted(glob.glob(os.path.join(self.bce_path, "CLUSTERS", "CL*", "cluster*.min.imaged.pdb"))) + clusters_order = self.order_clusters_min_distances(clusters) + dihedrals = self.list_all_dihedrals() + ordered_clusters = [clusters[x] for x in clusters_order] + for cluster in ordered_clusters: + self.calculate_cluster_angles(cluster, dihedrals, match_indexes=False) + + def _calculate_flexible_dihedrals(self): # here we rely on the output of the bce server to not change its # structrue in the future dihedral_file = glob.glob(os.path.join(self.bce_path, "DIHEDRALS", "*_new.ndx"))[0] @@ -47,7 +70,7 @@ def calculate(self): self.calculate_cluster_angles(cluster, dihedrals) - def calculate_cluster_angles(self, cluster_pdb, dihedral_list): + def calculate_cluster_angles(self, cluster_pdb, dihedral_list, match_indexes=True): """ Calculate dihedral angles from pdb @@ -55,9 +78,11 @@ def calculate_cluster_angles(self, cluster_pdb, dihedral_list): ---------- cluster_pdb: str Path to the cluster representative conformation - dihedral_list: list List of the tuples containing the atoms that form the dihedrals + match_indexes: bool + Whether to use the atom indices from the dihedral list or match to + the cluster structure before """ rdkit_wrapper = RDKitToolkitWrapper() pdb_dihedrals = [] @@ -66,10 +91,21 @@ def calculate_cluster_angles(self, cluster_pdb, dihedral_list): mol = molecule.Molecule(cluster_pdb, connectivity_template=self._molecule.rdkit_molecule) # we use substructure matching to ensure that the indices in the # clusters pdb and the input ligand are the same - rename_dict = {i: x for i, x in enumerate(rdkit_wrapper.get_substruct_match(self._molecule, mol))} + if match_indexes: + rename_dict = {i: x for i, x in enumerate(rdkit_wrapper.get_substruct_match(self._molecule, mol))} + else: + rename_dict = {i: x for i, x in enumerate(rdkit_wrapper.get_substruct_match(mol, self._molecule))} for dihedral in dihedral_list: - names = [self._topology.atoms[rename_dict[atom]].PDB_name for atom in dihedral] - angle = rdkit_wrapper.get_dihedral(mol, *dihedral, units="radians") + if match_indexes: + names = [self._topology.atoms[rename_dict[atom]].PDB_name for atom in dihedral] + angle = rdkit_wrapper.get_dihedral(mol, *dihedral, units="radians") + else: + names = [self._topology.atoms[atom].PDB_name for atom in dihedral] + if any(["H" in name for name in names]): + continue + dihedral_matched = tuple([rename_dict[index] for index in dihedral]) + angle = rdkit_wrapper.get_dihedral(mol, *dihedral_matched, units="radians") + pdb_dihedrals.append(names+[angle]) self.dihedral_library[cluster_pdb] = pdb_dihedrals @@ -93,6 +129,80 @@ def save(self, output_path): fw.write("ENDDIHEDRALS\n") fw.write("END") + def list_all_dihedrals(self): + dihedrals = [] + seen_dihedrals = set() + for proper in self._topology.propers: + dihedral_indexes = (proper.atom1_idx, proper.atom2_idx, proper.atom3_idx, proper.atom4_idx) + if dihedral_indexes in seen_dihedrals: + continue + seen_dihedrals.add(dihedral_indexes) + dihedrals.append(list(dihedral_indexes)) + return dihedrals + + + def order_clusters_min_distances(self, clusters): + """ + Order a list of structures in a way that jumps to the next structure + involves the minimum change + + Parameters + ---------- + + clusters: list + Path to the structures + """ + rdkit_wrapper = RDKitToolkitWrapper() + distances = np.zeros((len(clusters), len(clusters))) + cluster_molecules = [] + for cluster in clusters: + # use the input molecule as template since the cluster structures + # probably will not have proper stereochemistry + cluster_molecules.append(molecule.Molecule(cluster, connectivity_template=self._molecule.rdkit_molecule)) + for i, cluster_mol in enumerate(cluster_molecules): + for j, cluster_mol_2 in enumerate(cluster_molecules[i+1:], start=i+1): + rmsd_value = rdkit_wrapper.get_rmsd(cluster_mol, cluster_mol_2) + distances[i, j] = rmsd_value + distances[j, i] = rmsd_value + return find_optimal_path_from_matrix(distances) + + +def find_optimal_path_from_matrix(distances): + """ + Find a Minimum Spanning Tree from a distance matrix. It uses networkx to + create the graph and the MST + + Parameters + ---------- + + Returns + ------- + min_path : list + Nodes on the graph that minimise heuristically the total distance + """ + graph = nx.from_numpy_matrix(distances) + min_dist = 1e6 + min_path = None + for node in graph: + dist, path = find_heuristic_path(graph.copy(), distances, node) + if dist < min_dist: + min_dist = dist + min_path = path + return min_path + +def find_heuristic_path(graph, distances, start_node): + n = distances.shape[0] + path = [start_node] + dist = 0 + while len(path) < n: + node = path[-1] + new_neighbor = min(graph[node].items(), key=lambda x: x[1]["weight"]) + graph.remove_node(node) + path.append(new_neighbor[0]) + dist += new_neighbor[1]['weight'] + dist += distances[path[-1], start_node] + return dist, path + def read_dihedral_info_file(file_name): """ Read an ndx file containing the classification of all dihedrals and return @@ -102,6 +212,11 @@ def read_dihedral_info_file(file_name): ---------- file_name: str Path of the dihedrals file + + Returns + ------- + dihedrals : list + List of atom indices correponding to the flexible dihedrals """ read_dihedral = False dihedrals = [] diff --git a/peleffy/main.py b/peleffy/main.py index 7d7db030..63c41d71 100644 --- a/peleffy/main.py +++ b/peleffy/main.py @@ -58,6 +58,10 @@ def parse_args(args): parser.add_argument("--dihedrals_info_path", default=None, type=str, help="Path to the folder containing the BCE output" + " used to collect dihedral angles for PELE") + parser.add_argument("--dihedrals_mode", default="all_dihedrals", type=str, + help="Mode for extracting the dihedrals: 'all_dihedrals' will"+ + "extract the values fro all, while 'flexible_dihedrals' will extract"+ + "the values only for those that are marked as flexible") parser.add_argument('--with_solvent', dest='with_solvent', help="Generate solvent parameters for OBC", action='store_true') @@ -99,7 +103,7 @@ def run_peleffy(pdb_file, charge_method=DEFAULT_CHARGE_METHOD, exclude_terminal_rotamers=True, output=None, with_solvent=False, as_datalocal=False, - dihedral_path=None): + dihedral_path=None, dihedral_mode="all_dihedrals"): """ It runs peleffy. @@ -124,6 +128,10 @@ def run_peleffy(pdb_file, as_datalocal : bool Whether to save output files following PELE's DataLocal hierarchy or not + dihedral_path: str + Path to the BCE server outupt to use to extract dihedral angles + dihedral_mode: str + Select what kind of dihedrals to extract (all or only flexible) """ log = Logger() log.info('-' * 60) From 39e8c4049824e646fb627768d795f8d327d26951 Mon Sep 17 00:00:00 2001 From: jgilaber Date: Wed, 23 Dec 2020 22:25:40 +0100 Subject: [PATCH 16/27] Raise more informative error for missing BCEDihedrals --- peleffy/BCEDihedrals/BCEDihedrals.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/peleffy/BCEDihedrals/BCEDihedrals.py b/peleffy/BCEDihedrals/BCEDihedrals.py index e2223b13..0c62ee87 100644 --- a/peleffy/BCEDihedrals/BCEDihedrals.py +++ b/peleffy/BCEDihedrals/BCEDihedrals.py @@ -54,6 +54,8 @@ def calculate(self): def _calculate_all_dihedrals(self): clusters = sorted(glob.glob(os.path.join(self.bce_path, "CLUSTERS", "CL*", "cluster*.min.imaged.pdb"))) + if not clusters: + raise ValueError("Path to the BCE output does not contain a CLUSTERS folder, please check if the path is correct!") clusters_order = self.order_clusters_min_distances(clusters) dihedrals = self.list_all_dihedrals() ordered_clusters = [clusters[x] for x in clusters_order] From 8677dda4d1abe7cbc0328f6e6421b34749144053 Mon Sep 17 00:00:00 2001 From: jgilaber Date: Wed, 23 Dec 2020 22:26:50 +0100 Subject: [PATCH 17/27] Avoid calculating rotamer libraries with BCE --- peleffy/main.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/peleffy/main.py b/peleffy/main.py index 63c41d71..6d7866d2 100644 --- a/peleffy/main.py +++ b/peleffy/main.py @@ -172,8 +172,10 @@ def run_peleffy(pdb_file, output_path=output, as_datalocal=as_datalocal) - rotamer_library = peleffy.topology.RotamerLibrary(molecule) - rotamer_library.to_file(output_handler.get_rotamer_library_path()) + if dihedral_path is None: + # if dihedral_path is set, we don't want a rotamer library + rotamer_library = peleffy.topology.RotamerLibrary(molecule) + rotamer_library.to_file(output_handler.get_rotamer_library_path()) # Parameterize molecule with the selected force field log.info(' - Parameterizing molecule') @@ -208,6 +210,8 @@ def run_peleffy(pdb_file, log.info(' - All files were generated successfully:') log.info(' - {}'.format(output_handler.get_rotamer_library_path())) log.info(' - {}'.format(output_handler.get_impact_template_path())) + if dihedral_path is not None: + log.info(' - {}'.format(output_handler.get_dihedral_library_path())) if with_solvent: log.info(' - {}'.format(output_handler.get_solvent_template_path())) From 0bc819365c140b47f12574758f20b1d80aa7df65 Mon Sep 17 00:00:00 2001 From: jgilaber Date: Thu, 24 Dec 2020 11:01:57 +0100 Subject: [PATCH 18/27] Avoid logging rotamer library if not calculated --- peleffy/main.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/peleffy/main.py b/peleffy/main.py index 6d7866d2..894a5024 100644 --- a/peleffy/main.py +++ b/peleffy/main.py @@ -208,7 +208,8 @@ def run_peleffy(pdb_file, log.set_level(previous_level) log.info(' - All files were generated successfully:') - log.info(' - {}'.format(output_handler.get_rotamer_library_path())) + if dihedral_path is None: + log.info(' - {}'.format(output_handler.get_rotamer_library_path())) log.info(' - {}'.format(output_handler.get_impact_template_path())) if dihedral_path is not None: log.info(' - {}'.format(output_handler.get_dihedral_library_path())) From 06140f6eeb7e5eed9912e890638944292e5232cf Mon Sep 17 00:00:00 2001 From: jgilaber Date: Tue, 5 Jan 2021 14:11:13 +0100 Subject: [PATCH 19/27] Do not avoid reading dihedrals with hydrogens --- peleffy/BCEDihedrals/BCEDihedrals.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/peleffy/BCEDihedrals/BCEDihedrals.py b/peleffy/BCEDihedrals/BCEDihedrals.py index 0c62ee87..34110654 100644 --- a/peleffy/BCEDihedrals/BCEDihedrals.py +++ b/peleffy/BCEDihedrals/BCEDihedrals.py @@ -103,8 +103,6 @@ def calculate_cluster_angles(self, cluster_pdb, dihedral_list, match_indexes=Tru angle = rdkit_wrapper.get_dihedral(mol, *dihedral, units="radians") else: names = [self._topology.atoms[atom].PDB_name for atom in dihedral] - if any(["H" in name for name in names]): - continue dihedral_matched = tuple([rename_dict[index] for index in dihedral]) angle = rdkit_wrapper.get_dihedral(mol, *dihedral_matched, units="radians") From f0dfa8fbe2095971735f8ffe60be073549841d23 Mon Sep 17 00:00:00 2001 From: jgilaber Date: Wed, 17 Feb 2021 13:10:33 +0100 Subject: [PATCH 20/27] Add tests for outputHandler in the dihedral library --- peleffy/tests/test_utils.py | 49 +++++++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/peleffy/tests/test_utils.py b/peleffy/tests/test_utils.py index d0705a59..028e69f4 100644 --- a/peleffy/tests/test_utils.py +++ b/peleffy/tests/test_utils.py @@ -194,6 +194,9 @@ def test_non_datalocal_paths(self): assert output_handler.get_solvent_template_path() == \ './ligandParams.txt', \ 'Unexpected default solvent parameters path' + assert output_handler.get_dihedral_library_path() == \ + './{}.dihedral'.format(tag.upper()), \ + 'Unexpected default dihedral library path' # Initialize output handler with an output_path set with tempfile.TemporaryDirectory() as tmpdir: @@ -211,6 +214,9 @@ def test_non_datalocal_paths(self): assert output_handler.get_solvent_template_path() == \ os.path.join(tmpdir, 'output', 'ligandParams.txt'), \ 'Unexpected default solvent parameters path' + assert output_handler.get_dihedral_library_path() == \ + os.path.join(tmpdir, 'output', '{}.dihedral'.format(tag.upper())), \ + 'Unexpected default dihedral library path' def test_datalocal_paths_for_openff(self): """It tests the datalocal paths assignment for OpenFF.""" @@ -244,6 +250,10 @@ def test_datalocal_paths_for_openff(self): create_missing_folders=False) == \ './DataLocal/OBC/ligandParams.txt', \ 'Unexpected default solvent parameters path' + assert output_handler.get_dihedral_library_path( + create_missing_folders=False) == \ + './DataLocal/Dihedrals/'+'{}.dihedral'.format(tag.upper()), \ + 'Unexpected default dihedral library path' with tempfile.TemporaryDirectory() as tmpdir: output_handler = OutputPathHandler( @@ -265,6 +275,10 @@ def test_datalocal_paths_for_openff(self): create_missing_folders=False) == \ tmpdir + '/output/DataLocal/OBC/ligandParams.txt', \ 'Unexpected default solvent parameters path' + assert output_handler.get_dihedral_library_path( + create_missing_folders=False) == \ + tmpdir + '/output/DataLocal/Dihedrals/'+'{}.dihedral'.format(tag.upper()), \ + 'Unexpected default dihedral library path' def test_datalocal_paths_for_opls(self): """It tests the datalocal paths assignment for OPLS2005.""" @@ -299,6 +313,10 @@ def test_datalocal_paths_for_opls(self): create_missing_folders=False) == \ './DataLocal/OBC/ligandParams.txt', \ 'Unexpected default solvent parameters path' + assert output_handler.get_dihedral_library_path( + create_missing_folders=False) == \ + './DataLocal/Dihedrals/'+'{}.dihedral'.format(tag.upper()), \ + 'Unexpected default dihedral library path' # Initialize output handler with an output_path set with tempfile.TemporaryDirectory() as tmpdir: @@ -321,6 +339,10 @@ def test_datalocal_paths_for_opls(self): os.path.join(tmpdir, 'output', 'DataLocal/OBC/ligandParams.txt'), \ 'Unexpected default solvent parameters path' + assert output_handler.get_dihedral_library_path( + create_missing_folders=False) == \ + os.path.join(tmpdir, 'output/DataLocal/Dihedrals/'+'{}.dihedral'.format(tag.upper())), \ + 'Unexpected default dihedral library path' def test_datalocal_paths_for_offopls(self): """ @@ -360,6 +382,10 @@ def test_datalocal_paths_for_offopls(self): create_missing_folders=False) == \ './DataLocal/OBC/ligandParams.txt', \ 'Unexpected default solvent parameters path' + assert output_handler.get_dihedral_library_path( + create_missing_folders=False) == \ + './DataLocal/Dihedrals/'+'{}.dihedral'.format(tag.upper()), \ + 'Unexpected default dihedral library path' # Initialize output handler with an output_path set with tempfile.TemporaryDirectory() as tmpdir: @@ -382,6 +408,10 @@ def test_datalocal_paths_for_offopls(self): os.path.join(tmpdir, 'output', 'DataLocal/OBC/ligandParams.txt'), \ 'Unexpected default solvent parameters path' + assert output_handler.get_dihedral_library_path( + create_missing_folders=False) == \ + os.path.join(tmpdir, 'output/DataLocal/Dihedrals/'+'{}.dihedral'.format(tag.upper())), \ + 'Unexpected default dihedral library path' # Initialize output handler without output_path output_handler = OutputPathHandler(molecule, hybridff, @@ -405,6 +435,10 @@ def test_datalocal_paths_for_offopls(self): create_missing_folders=False) == \ './DataLocal/OBC/ligandParams.txt', \ 'Unexpected default solvent parameters path' + assert output_handler.get_dihedral_library_path( + create_missing_folders=False) == \ + './DataLocal/Dihedrals/'+'{}.dihedral'.format(tag.upper()), \ + 'Unexpected default dihedral library path' # Initialize output handler with an output_path set with tempfile.TemporaryDirectory() as tmpdir: @@ -427,6 +461,11 @@ def test_datalocal_paths_for_offopls(self): os.path.join(tmpdir, 'output', 'DataLocal/OBC/ligandParams.txt'), \ 'Unexpected default solvent parameters path' + assert output_handler.get_dihedral_library_path( + create_missing_folders=False) == \ + os.path.join(tmpdir, 'output', + 'DataLocal/Dihedrals/'+'{}.dihedral'.format(tag.upper())), \ + 'Unexpected default dihedral library path' def test_folder_creation(self): """ @@ -464,6 +503,11 @@ def test_folder_creation(self): path_dir = os.path.dirname(path) assert os.path.isdir(path_dir) is False, \ 'This directory should not exist' + path = output_handler.get_dihedral_library_path( + create_missing_folders=False) + path_dir = os.path.dirname(path) + assert os.path.isdir(path_dir) is False, \ + 'This directory should not exist' # Initialize output handler with an output_path set with tempfile.TemporaryDirectory() as tmpdir: @@ -489,6 +533,11 @@ def test_folder_creation(self): path_dir = os.path.dirname(path) assert os.path.isdir(path_dir) is True, \ 'This directory should exist' + path = output_handler.get_dihedral_library_path( + create_missing_folders=True) + path_dir = os.path.dirname(path) + assert os.path.isdir(path_dir) is True, \ + 'This directory should exist' class TestMAEParser(object): """ From 8d1b214aacf5ec80e755ed3663ed699744dad873 Mon Sep 17 00:00:00 2001 From: jgilaber Date: Wed, 17 Feb 2021 14:57:30 +0100 Subject: [PATCH 21/27] Pass dihedrals_mode argument to main --- peleffy/main.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/peleffy/main.py b/peleffy/main.py index 6a9c474a..95bc7ed8 100644 --- a/peleffy/main.py +++ b/peleffy/main.py @@ -60,7 +60,7 @@ def parse_args(args): + " used to collect dihedral angles for PELE") parser.add_argument("--dihedrals_mode", default="all_dihedrals", type=str, help="Mode for extracting the dihedrals: 'all_dihedrals' will"+ - "extract the values fro all, while 'flexible_dihedrals' will extract"+ + "extract the values for all, while 'flexible_dihedrals' will extract"+ "the values only for those that are marked as flexible") parser.add_argument('--with_solvent', dest='with_solvent', help="Generate solvent parameters for OBC", @@ -281,7 +281,8 @@ def main(args): with_solvent=args.with_solvent, as_datalocal=args.as_datalocal, dihedral_path=args.dihedrals_info_path, - charges_from_file=args.charges_from_file) + charges_from_file=args.charges_from_file, + dihedral_mode=args.dihedrals_mode) if __name__ == '__main__': From 96624309c8b6240ee1ae97db12ccac1d1a29ccbf Mon Sep 17 00:00:00 2001 From: jgilaber Date: Wed, 17 Feb 2021 14:57:47 +0100 Subject: [PATCH 22/27] Test parsing dihedral library parameters in main --- peleffy/tests/test_main.py | 62 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) diff --git a/peleffy/tests/test_main.py b/peleffy/tests/test_main.py index b6e2e57f..cd4b3ae8 100644 --- a/peleffy/tests/test_main.py +++ b/peleffy/tests/test_main.py @@ -81,6 +81,10 @@ def test_peleffy_argparse(self): 'Unexpected silent settings were parsed' assert parsed_args.with_solvent is False, \ 'Unexpected with_solvent settings were parsed' + assert parsed_args.dihedrals_info_path is None, \ + 'Unexpected dihedral_path settings were passed' + assert parsed_args.dihedrals_mode == "all_dihedrals", \ + 'Unexpected dihedrals_mode settings were passed' # Test custom shorts parsed_args = parse_args(['toluene.pdb', @@ -109,6 +113,10 @@ def test_peleffy_argparse(self): 'Unexpected silent settings were parsed' assert parsed_args.with_solvent is False, \ 'Unexpected with_solvent settings were parsed' + assert parsed_args.dihedrals_info_path is None, \ + 'Unexpected dihedral_path settings were passed' + assert parsed_args.dihedrals_mode == "all_dihedrals", \ + 'Unexpected dihedrals_mode settings were passed' # Test custom longs parsed_args = parse_args(['methane.pdb', @@ -138,6 +146,10 @@ def test_peleffy_argparse(self): 'Unexpected silent settings were parsed' assert parsed_args.with_solvent is False, \ 'Unexpected with_solvent settings were parsed' + assert parsed_args.dihedrals_info_path is None, \ + 'Unexpected dihedral_path settings were passed' + assert parsed_args.dihedrals_mode == "all_dihedrals", \ + 'Unexpected dihedrals_mode settings were passed' # Test unexpected charge method with pytest.raises(SystemExit) as pytest_wrapped_e: @@ -170,6 +182,10 @@ def test_peleffy_argparse(self): 'Unexpected silent settings were parsed' assert parsed_args.with_solvent is False, \ 'Unexpected with_solvent settings were parsed' + assert parsed_args.dihedrals_info_path is None, \ + 'Unexpected dihedral_path settings were passed' + assert parsed_args.dihedrals_mode == "all_dihedrals", \ + 'Unexpected dihedrals_mode settings were passed' # Test charges_from_file argument parsed_args = parse_args(['BHP.pdb', @@ -197,6 +213,10 @@ def test_peleffy_argparse(self): 'Unexpected silent settings were parsed' assert parsed_args.with_solvent is False, \ 'Unexpected with_solvent settings were parsed' + assert parsed_args.dihedrals_info_path is None, \ + 'Unexpected dihedral_path settings were passed' + assert parsed_args.dihedrals_mode == "all_dihedrals", \ + 'Unexpected dihedrals_mode settings were passed' # Test include_terminal_rotamers argument parsed_args = parse_args(['methane.pdb', @@ -222,6 +242,10 @@ def test_peleffy_argparse(self): 'Unexpected silent settings were parsed' assert parsed_args.with_solvent is False, \ 'Unexpected with_solvent settings were parsed' + assert parsed_args.dihedrals_info_path is None, \ + 'Unexpected dihedral_path settings were passed' + assert parsed_args.dihedrals_mode == "all_dihedrals", \ + 'Unexpected dihedrals_mode settings were passed' # Test silent argument parsed_args = parse_args(['methane.pdb', @@ -247,6 +271,10 @@ def test_peleffy_argparse(self): 'Unexpected silent settings were parsed' assert parsed_args.with_solvent is False, \ 'Unexpected with_solvent settings were parsed' + assert parsed_args.dihedrals_info_path is None, \ + 'Unexpected dihedral_path settings were passed' + assert parsed_args.dihedrals_mode == "all_dihedrals", \ + 'Unexpected dihedrals_mode settings were passed' parse_args(['methane.pdb', '-s']) == parse_args(['methane.pdb', '--silent']) @@ -275,10 +303,44 @@ def test_peleffy_argparse(self): 'Unexpected silent settings were parsed' assert parsed_args.with_solvent is False, \ 'Unexpected with_solvent settings were parsed' + assert parsed_args.dihedrals_info_path is None, \ + 'Unexpected dihedral_path settings were passed' + assert parsed_args.dihedrals_mode == "all_dihedrals", \ + 'Unexpected dihedrals_mode settings were passed' parse_args(['methane.pdb', '-d']) == parse_args(['methane.pdb', '--debug']) + # Test dihedral library arguments + parsed_args = parse_args(['methane.pdb', + '--dihedrals_info_path', 'test_path', + '--dihedrals_mode', 'flexible_dihedrals']) + + assert parsed_args.as_datalocal is False, \ + 'Unexpected as_datalocal settings were parsed' + assert parsed_args.charge_method == 'am1bcc', \ + 'Unexpected charge_method settings were parsed' + assert parsed_args.debug is False, \ + 'Unexpected debug settings were parsed' + assert parsed_args.forcefield == 'openff_unconstrained-1.3.0.offxml', \ + 'Unexpected forcefield settings were parsed' + assert parsed_args.include_terminal_rotamers is False, \ + 'Unexpected include_terminal_rotamers settings were parsed' + assert parsed_args.output is None, \ + 'Unexpected output settings were parsed' + assert parsed_args.pdb_file == 'methane.pdb', \ + 'Unexpected pdb_file settings were parsed' + assert parsed_args.resolution == 30, \ + 'Unexpected resolution settings were parsed' + assert parsed_args.silent is False, \ + 'Unexpected silent settings were parsed' + assert parsed_args.with_solvent is False, \ + 'Unexpected with_solvent settings were parsed' + assert parsed_args.dihedrals_info_path == "test_path", \ + 'Unexpected dihedral_path settings were passed' + assert parsed_args.dihedrals_mode == "flexible_dihedrals", \ + 'Unexpected dihedrals_mode settings were passed' + def test_peleffy_main(self): """It checks the main function of peleffy.""" from peleffy.main import parse_args, main From 88754d395871caa7d9eb0dfc922829e6ddd62f02 Mon Sep 17 00:00:00 2001 From: jgilaber Date: Wed, 17 Feb 2021 15:50:16 +0100 Subject: [PATCH 23/27] Add tests for new functions in rdkit toolkit --- .../data/ligands/trimethylglycine_moved.pdb | 47 ++++++++++++++++++ peleffy/tests/test_toolkits.py | 48 +++++++++++++++++++ 2 files changed, 95 insertions(+) create mode 100644 peleffy/data/ligands/trimethylglycine_moved.pdb diff --git a/peleffy/data/ligands/trimethylglycine_moved.pdb b/peleffy/data/ligands/trimethylglycine_moved.pdb new file mode 100644 index 00000000..bc725155 --- /dev/null +++ b/peleffy/data/ligands/trimethylglycine_moved.pdb @@ -0,0 +1,47 @@ +REMARK 4 COMPLIES WITH FORMAT V. 3.0, 1-DEC-2006 +REMARK 888 +REMARK 888 WRITTEN BY MAESTRO (A PRODUCT OF SCHRODINGER, LLC) +TITLE trimethylglycine_2 +MODEL 1 +HETATM 1 C1 UNL 1 0.028 -1.392 0.008 1.00 0.00 C +HETATM 2 N1 UNL 1 -0.273 0.003 -0.010 1.00 0.00 N1+ +HETATM 3 C2 UNL 1 -0.985 0.255 -1.275 1.00 0.00 C +HETATM 4 C3 UNL 1 -1.213 0.331 1.044 1.00 0.00 C +HETATM 5 C4 UNL 1 0.855 0.875 -0.041 1.00 0.00 C +HETATM 6 C5 UNL 1 2.061 0.381 0.593 1.00 0.00 C +HETATM 7 O1 UNL 1 3.096 1.084 0.610 1.00 0.00 O +HETATM 8 O2 UNL 1 2.145 -0.855 1.201 1.00 0.00 O1- +HETATM 9 H1 UNL 1 0.016 -1.836 -1.009 1.00 0.00 H +HETATM 10 H2 UNL 1 -0.799 -1.910 0.570 1.00 0.00 H +HETATM 11 H3 UNL 1 0.958 -1.673 0.486 1.00 0.00 H +HETATM 12 H4 UNL 1 -0.224 0.185 -2.088 1.00 0.00 H +HETATM 13 H5 UNL 1 -1.834 -0.417 -1.403 1.00 0.00 H +HETATM 14 H6 UNL 1 -1.314 1.314 -1.267 1.00 0.00 H +HETATM 15 H7 UNL 1 -2.179 0.700 0.673 1.00 0.00 H +HETATM 16 H8 UNL 1 -1.429 -0.586 1.636 1.00 0.00 H +HETATM 17 H9 UNL 1 -0.809 1.096 1.743 1.00 0.00 H +HETATM 18 H10 UNL 1 1.028 1.180 -1.092 1.00 0.00 H +HETATM 19 H11 UNL 1 0.525 1.823 0.481 1.00 0.00 H +CONECT 1 2 9 10 11 +CONECT 2 1 3 4 5 +CONECT 3 2 12 13 14 +CONECT 4 2 15 16 17 +CONECT 5 2 6 18 19 +CONECT 6 5 7 8 +CONECT 6 7 +CONECT 7 6 +CONECT 7 6 +CONECT 8 6 +CONECT 9 1 +CONECT 10 1 +CONECT 11 1 +CONECT 12 3 +CONECT 13 3 +CONECT 14 3 +CONECT 15 4 +CONECT 16 4 +CONECT 17 4 +CONECT 18 5 +CONECT 19 5 +ENDMDL +END diff --git a/peleffy/tests/test_toolkits.py b/peleffy/tests/test_toolkits.py index 2ec31ad6..1e460131 100644 --- a/peleffy/tests/test_toolkits.py +++ b/peleffy/tests/test_toolkits.py @@ -2,6 +2,7 @@ This module contains the tests to check the peleffy's toolkits. """ +import numpy as np import pytest from peleffy.forcefield.parameters import OPLS2005ParameterWrapper @@ -270,3 +271,50 @@ def test_pdb_parsers(self): block2 = MolToPDBBlock(rdkit_mol2) assert block1 == block2, 'Unexpected pair of RDKit molecules' + + def test_dihedral_angle(self): + """It checks that the dihedral angle calculator works well.""" + + from peleffy.topology import Molecule + from peleffy.utils.toolkits import RDKitToolkitWrapper + from peleffy.utils import get_data_file_path + + wrapper = RDKitToolkitWrapper() + + pdb_path = get_data_file_path('ligands/trimethylglycine.pdb') + m = Molecule(pdb_path) + dihedral_degrees = wrapper.get_dihedral(m, 2, 1, 4, 5, units="degrees") + dihedral_rad = wrapper.get_dihedral(m, 2, 1, 4, 5) + np.testing.assert_almost_equal(dihedral_degrees, -176.348, decimal=2) + np.testing.assert_almost_equal(dihedral_degrees, np.rad2deg(dihedral_rad), decimal=3) + + def test_dihedral_angle_2(self): + """It checks that the dihedral angle calculator works well.""" + + from peleffy.topology import Molecule + from peleffy.utils.toolkits import RDKitToolkitWrapper + from peleffy.utils import get_data_file_path + + wrapper = RDKitToolkitWrapper() + + pdb_path = get_data_file_path('ligands/trimethylglycine.pdb') + m = Molecule(pdb_path) + dihedral_degrees = wrapper.get_dihedral(m, 17, 4, 5, 6, units="degrees") + dihedral_rad = wrapper.get_dihedral(m, 17, 4, 5, 6) + np.testing.assert_almost_equal(dihedral_degrees, 54.828, decimal=2) + np.testing.assert_almost_equal(dihedral_degrees, np.rad2deg(dihedral_rad), decimal=3) + + def test_rmsd(self): + """It checks that the rmsd calculator works well.""" + + from peleffy.topology import Molecule + from peleffy.utils.toolkits import RDKitToolkitWrapper + from peleffy.utils import get_data_file_path + + wrapper = RDKitToolkitWrapper() + + pdb_path = get_data_file_path('ligands/trimethylglycine.pdb') + m = Molecule(pdb_path) + pdb_path2 = get_data_file_path('ligands/trimethylglycine_moved.pdb') + m2 = Molecule(pdb_path2) + np.testing.assert_almost_equal(wrapper.get_rmsd(m, m2), 0.3346, decimal=3) From d118a9beaf29eab703688ff3282886f2e982ce1f Mon Sep 17 00:00:00 2001 From: jgilaber Date: Wed, 17 Feb 2021 17:31:11 +0100 Subject: [PATCH 24/27] Add tests for BCEDihedrals module --- peleffy/BCEDihedrals/BCEDihedrals.py | 2 +- peleffy/data/parameters/ETH.dihedral | 10 ++++ peleffy/tests/test_BCEDihedrals.py | 71 ++++++++++++++++++++++++++++ peleffy/tests/test_integrity.py | 1 + 4 files changed, 83 insertions(+), 1 deletion(-) create mode 100644 peleffy/data/parameters/ETH.dihedral create mode 100644 peleffy/tests/test_BCEDihedrals.py diff --git a/peleffy/BCEDihedrals/BCEDihedrals.py b/peleffy/BCEDihedrals/BCEDihedrals.py index 34110654..3fc454cf 100644 --- a/peleffy/BCEDihedrals/BCEDihedrals.py +++ b/peleffy/BCEDihedrals/BCEDihedrals.py @@ -127,7 +127,7 @@ def save(self, output_path): for dihedral in dihedral_values: fw.write("{:s} {:s} {:s} {:s} {:5f}\n".format(*dihedral)) fw.write("ENDDIHEDRALS\n") - fw.write("END") + fw.write("END\n") def list_all_dihedrals(self): dihedrals = [] diff --git a/peleffy/data/parameters/ETH.dihedral b/peleffy/data/parameters/ETH.dihedral new file mode 100644 index 00000000..f5979903 --- /dev/null +++ b/peleffy/data/parameters/ETH.dihedral @@ -0,0 +1,10 @@ +* DIHEDRAL LIBRARY FILE +* File generated with peleffy-1.0.0+164.g06140f6 +* File: ../../PELEPipelineSanja/4S0V-LIG/CLUSTERS/CL1/cluster1.min.imaged.pdb +UNK 4 1 +_H1_ _C1_ _C2_ _H3_ 0.001146 +_H1_ _C1_ _C2_ _H4_ -3.141280 +_H2_ _C1_ _C2_ _H3_ -3.141120 +_H2_ _C1_ _C2_ _H4_ -0.000361 +ENDDIHEDRALS +END diff --git a/peleffy/tests/test_BCEDihedrals.py b/peleffy/tests/test_BCEDihedrals.py new file mode 100644 index 00000000..c37052ed --- /dev/null +++ b/peleffy/tests/test_BCEDihedrals.py @@ -0,0 +1,71 @@ +""" +This module contains tests that check that the BCEDihedrals module. +""" + +import numpy as np +import pytest + +def build_mock_BCEDihedrals(pdb_file): + from peleffy.topology import Molecule + from peleffy.BCEDihedrals import BCEDihedrals + from peleffy.forcefield import ForceFieldSelector + from peleffy.topology import Topology + + molecule = Molecule(pdb_file) + ff_selector = ForceFieldSelector() + forcefield = ff_selector.get_by_name("opls2005") + parameters = forcefield.parameterize(molecule, charge_method="opls2005") + topology = Topology(molecule, parameters) + return BCEDihedrals(topology, "", mode="") + +class TestBCEDihderals(object): + """BCEDihedrals test.""" + + from peleffy.BCEDihedrals import BCEDihedrals + + def test_list_dihedrals(self): + from peleffy.utils import get_data_file_path + + pdb_path = get_data_file_path('ligands/ethylene.pdb') + golden_dihedrals = [(2, 0, 1, 4), (2, 0, 1, 5), (3, 0, 1, 4), (3, 0, 1, 5)] + bce_obj = build_mock_BCEDihedrals(pdb_path) + dihedrals_indices = bce_obj.list_all_dihedrals() + dihedrals_indices.sort() + np.testing.assert_array_equal(dihedrals_indices, golden_dihedrals) + + def test_calculate_cluster_angles(self): + from peleffy.utils import get_data_file_path + + golden_dihedrals = [(2, 0, 1, 4), (2, 0, 1, 5), (3, 0, 1, 4), (3, 0, 1, 5)] + pdb_path = get_data_file_path('ligands/ethylene.pdb') + golden_angles = [ ["_H1_", "_C1_", "_C2_", "_H3_", 0.001151], + ["_H1_", "_C1_", "_C2_", "_H4_", -3.141278], + ["_H2_", "_C1_", "_C2_", "_H3_", -3.141278], + ["_H2_", "_C1_", "_C2_", "_H4_", -0.000366]] + bce_obj = build_mock_BCEDihedrals(pdb_path) + bce_obj.calculate_cluster_angles(pdb_path, golden_dihedrals, match_indexes=False) + for dih1, dih2 in zip(bce_obj.dihedral_library[pdb_path], golden_angles): + np.testing.assert_array_equal(dih1[:4], dih2[:4]) + np.testing.assert_almost_equal(dih1[-1], dih2[-1], decimal=3) + + def test_write_dihedral_library(self): + import os + import tempfile + from peleffy.utils import get_data_file_path, temporary_cd + + golden_dihedrals = [(2, 0, 1, 4), (2, 0, 1, 5), (3, 0, 1, 4), (3, 0, 1, 5)] + pdb_path = get_data_file_path('ligands/ethylene.pdb') + bce_obj = build_mock_BCEDihedrals(pdb_path) + bce_obj.calculate_cluster_angles(pdb_path, golden_dihedrals, match_indexes=False) + calculated_lines = [] + with tempfile.TemporaryDirectory() as tmpdir: + with temporary_cd(tmpdir): + bce_obj.save(os.path.join(tmpdir, "ETH.dihedral")) + with open(os.path.join(tmpdir, "ETH.dihedral")) as f: + calculated_lines = f.readlines() + calculated_lines = calculated_lines[3:] + golden_dihedral_library_path = get_data_file_path('parameters/ETH.dihedral') + with open(golden_dihedral_library_path) as f: + golden_lines = f.readlines()[3:] + + assert golden_lines == calculated_lines diff --git a/peleffy/tests/test_integrity.py b/peleffy/tests/test_integrity.py index 475dad89..f1871fd4 100644 --- a/peleffy/tests/test_integrity.py +++ b/peleffy/tests/test_integrity.py @@ -17,6 +17,7 @@ def test_modules(self): from peleffy import template from peleffy import topology from peleffy import utils + from peleffy import BCEDihedrals except ImportError as e: raise AssertionError('The following peleffy module is missing: ' From 68e206cbc9782870e10b7fceffe44f0146a4070a Mon Sep 17 00:00:00 2001 From: jgilaber Date: Thu, 11 Mar 2021 17:01:43 +0100 Subject: [PATCH 25/27] Switch from calculating dihedrals to translations In this new mode, the tranlation from each atom to the root atom of the template is calculated instead of the dihedral angle --- peleffy/BCEConformations/BCEConformations.py | 170 +++++++++++++ peleffy/BCEConformations/__init__.py | 1 + peleffy/BCEDihedrals/BCEDihedrals.py | 237 ------------------- peleffy/BCEDihedrals/__init__.py | 1 - peleffy/main.py | 37 ++- peleffy/utils/output.py | 20 +- 6 files changed, 197 insertions(+), 269 deletions(-) create mode 100644 peleffy/BCEConformations/BCEConformations.py create mode 100644 peleffy/BCEConformations/__init__.py delete mode 100644 peleffy/BCEDihedrals/BCEDihedrals.py delete mode 100644 peleffy/BCEDihedrals/__init__.py diff --git a/peleffy/BCEConformations/BCEConformations.py b/peleffy/BCEConformations/BCEConformations.py new file mode 100644 index 00000000..391e8f61 --- /dev/null +++ b/peleffy/BCEConformations/BCEConformations.py @@ -0,0 +1,170 @@ +""" +This module contains classes and functions involved in the +creation of a library of dihedral angles from the output of the +BCE server +""" +import os +import glob +import numpy as np +import networkx as nx +import peleffy +from peleffy.utils import Logger +from peleffy.template import impact +from peleffy.topology import molecule +from peleffy.utils.toolkits import RDKitToolkitWrapper + +class BCEConformations(object): + """ + A class to produce a library of conformations from the output + of the BCE server + """ + def __init__(self, topology_obj, bce_output_path, verbose=False): + """ + Initializes a BCEConformations object + + Parameters + ---------- + topology_obj : An peleffy.topology.Topology + A Topology object that contains the ligand's information + bce_output_path: str + Path where the output from the BCE server is stored + + """ + self._topology = topology_obj + self._molecule = topology_obj.molecule + self.bce_path = bce_output_path + self.conformations_library = {} + self.verbose = verbose + + def calculate(self): + """ + Calculate conformation library from the bce output + """ + logger = Logger() + logger.info(' - Calculating conformation library') + if not self.verbose: + logger.set_level('WARNING') + self._calculate_all_conformations() + + def _calculate_all_conformations(self): + clusters = sorted(glob.glob(os.path.join(self.bce_path, "CLUSTERS", "CL*", "cluster*.min.imaged.pdb"))) + if not clusters: + raise ValueError("Path to the BCE output does not contain a CLUSTERS folder, please check if the path is correct!") + clusters_order = self.order_clusters_min_distances(clusters) + ordered_clusters = [clusters[x] for x in clusters_order] + for cluster in ordered_clusters: + self.calculate_cluster_offsets(cluster) + + + def calculate_cluster_offsets(self, cluster_pdb): + """ + Calculate dihedral angles from pdb + + Parameters + ---------- + cluster_pdb: str + Path to the cluster representative conformation + """ + rdkit_wrapper = RDKitToolkitWrapper() + conformation_offsets = [] + # use the input molecule as template since the cluster structures + # probably will not have proper stereochemistry + mol = molecule.Molecule(cluster_pdb, connectivity_template=self._molecule.rdkit_molecule) + cluster_coordinates = mol.get_conformer() + topology_to_cluster = {i: x for i, x in enumerate(rdkit_wrapper.get_substruct_match(mol, self._molecule))} + cluster_to_topology = {v: k for k, v in topology_to_cluster.items()} + template = impact.Impact(self._topology) + root_coordinates = cluster_coordinates[topology_to_cluster[find_index_root(template.topology, self._topology)]] + for i, coords in enumerate(cluster_coordinates): + conformation_offsets.append([self._topology.atoms[cluster_to_topology[i]].PDB_name]+(coords-root_coordinates).tolist()) + self.conformations_library[cluster_pdb] = conformation_offsets + + def save(self, output_path): + """ + Save the conformation library + + Parameters + ---------- + output_path: str + Where to save the conformation library + """ + with open(output_path, "w") as fw: + fw.write("* CONFORMATIONS LIBRARY FILE\n") + fw.write('* File generated with peleffy-{}\n'.format( peleffy.__version__)) + for conformation_file, conformation_values in self.conformations_library.items(): + fw.write("* File: {:s}\n".format(conformation_file)) + fw.write("{:s} {:d} {:d}\n".format(self._molecule.tag.upper(), len(conformation_values), len(self.conformations_library))) + for values in conformation_values: + fw.write("{:s} {:5f} {:5f} {:5f}\n".format(*values)) + fw.write("ENDCONFORMATION\n") + fw.write("END\n") + + def order_clusters_min_distances(self, clusters): + """ + Order a list of structures in a way that jumps to the next structure + involves the minimum change + + Parameters + ---------- + + clusters: list + Path to the structures + """ + rdkit_wrapper = RDKitToolkitWrapper() + distances = np.zeros((len(clusters), len(clusters))) + cluster_molecules = [] + for cluster in clusters: + # use the input molecule as template since the cluster structures + # probably will not have proper stereochemistry + cluster_molecules.append(molecule.Molecule(cluster, connectivity_template=self._molecule.rdkit_molecule)) + for i, cluster_mol in enumerate(cluster_molecules): + for j, cluster_mol_2 in enumerate(cluster_molecules[i+1:], start=i+1): + rmsd_value = rdkit_wrapper.get_rmsd(cluster_mol, cluster_mol_2) + distances[i, j] = rmsd_value + distances[j, i] = rmsd_value + return find_optimal_path_from_matrix(distances) + + +def find_optimal_path_from_matrix(distances): + """ + Find a Minimum Spanning Tree from a distance matrix. It uses networkx to + create the graph and the MST + + Parameters + ---------- + + Returns + ------- + min_path : list + Nodes on the graph that minimise heuristically the total distance + """ + graph = nx.from_numpy_matrix(distances) + min_dist = 1e6 + min_path = None + for node in graph: + dist, path = find_heuristic_path(graph.copy(), distances, node) + if dist < min_dist: + min_dist = dist + min_path = path + return min_path + +def find_heuristic_path(graph, distances, start_node): + n = distances.shape[0] + path = [start_node] + dist = 0 + while len(path) < n: + node = path[-1] + new_neighbor = min(graph[node].items(), key=lambda x: x[1]["weight"]) + graph.remove_node(node) + path.append(new_neighbor[0]) + dist += new_neighbor[1]['weight'] + dist += distances[path[-1], start_node] + return dist, path + + +def find_index_root(sorted_topology, topology): + atom_name = sorted_topology.atoms[0].PDB_name + for i, atom in enumerate(topology.atoms): + if atom.PDB_name == atom_name: + return i + return None diff --git a/peleffy/BCEConformations/__init__.py b/peleffy/BCEConformations/__init__.py new file mode 100644 index 00000000..ad12d12c --- /dev/null +++ b/peleffy/BCEConformations/__init__.py @@ -0,0 +1 @@ +from .BCEConformations import BCEConformations diff --git a/peleffy/BCEDihedrals/BCEDihedrals.py b/peleffy/BCEDihedrals/BCEDihedrals.py deleted file mode 100644 index 3fc454cf..00000000 --- a/peleffy/BCEDihedrals/BCEDihedrals.py +++ /dev/null @@ -1,237 +0,0 @@ -""" -This module contains classes and functions involved in the -creation of a library of dihedral angles from the output of the -BCE server -""" -import os -import glob -import numpy as np -import networkx as nx -import peleffy -from peleffy.utils import Logger -from peleffy.topology import molecule -from peleffy.utils.toolkits import RDKitToolkitWrapper - -class BCEDihedrals(object): - """ - A class to produce a library of dihedral angles from the output - of the BCE server - """ - def __init__(self, topology_obj, bce_output_path, mode="all_dihedrals", verbose=False): - """ - Initializes a BCEDihedrals object - - Parameters - ---------- - topology_obj : An peleffy.topology.Topology - A Topology object that contains the ligand's information - bce_output_path: str - Path where the output from the BCE server is stored - mode: str - Whether to extract all dihedrals or only those marked as flexible - - """ - self._topology = topology_obj - self._molecule = topology_obj.molecule - self.bce_path = bce_output_path - self.dihedral_library = {} - self.mode = mode - self.verbose = verbose - - def calculate(self): - """ - Calculate dihedrals library from the bce output - """ - logger = Logger() - logger.info(' - Calculating dihedral library') - logger.info(' - with mode: {}'.format(self.mode)) - if not self.verbose: - logger.set_level('WARNING') - if self.mode == "all_dihedrals": - self._calculate_all_dihedrals() - elif self.mode == "flexible_dihedrals": - self._calculate_flexible_dihedrals() - - def _calculate_all_dihedrals(self): - clusters = sorted(glob.glob(os.path.join(self.bce_path, "CLUSTERS", "CL*", "cluster*.min.imaged.pdb"))) - if not clusters: - raise ValueError("Path to the BCE output does not contain a CLUSTERS folder, please check if the path is correct!") - clusters_order = self.order_clusters_min_distances(clusters) - dihedrals = self.list_all_dihedrals() - ordered_clusters = [clusters[x] for x in clusters_order] - for cluster in ordered_clusters: - self.calculate_cluster_angles(cluster, dihedrals, match_indexes=False) - - def _calculate_flexible_dihedrals(self): - # here we rely on the output of the bce server to not change its - # structrue in the future - dihedral_file = glob.glob(os.path.join(self.bce_path, "DIHEDRALS", "*_new.ndx"))[0] - dihedrals = read_dihedral_info_file(dihedral_file) - clusters = sorted(glob.glob(os.path.join(self.bce_path, "CLUSTERS", "CL*", "cluster*.min.imaged.pdb"))) - for cluster in clusters: - self.calculate_cluster_angles(cluster, dihedrals) - - - def calculate_cluster_angles(self, cluster_pdb, dihedral_list, match_indexes=True): - """ - Calculate dihedral angles from pdb - - Parameters - ---------- - cluster_pdb: str - Path to the cluster representative conformation - dihedral_list: list - List of the tuples containing the atoms that form the dihedrals - match_indexes: bool - Whether to use the atom indices from the dihedral list or match to - the cluster structure before - """ - rdkit_wrapper = RDKitToolkitWrapper() - pdb_dihedrals = [] - # use the input molecule as template since the cluster structures - # probably will not have proper stereochemistry - mol = molecule.Molecule(cluster_pdb, connectivity_template=self._molecule.rdkit_molecule) - # we use substructure matching to ensure that the indices in the - # clusters pdb and the input ligand are the same - if match_indexes: - rename_dict = {i: x for i, x in enumerate(rdkit_wrapper.get_substruct_match(self._molecule, mol))} - else: - rename_dict = {i: x for i, x in enumerate(rdkit_wrapper.get_substruct_match(mol, self._molecule))} - for dihedral in dihedral_list: - if match_indexes: - names = [self._topology.atoms[rename_dict[atom]].PDB_name for atom in dihedral] - angle = rdkit_wrapper.get_dihedral(mol, *dihedral, units="radians") - else: - names = [self._topology.atoms[atom].PDB_name for atom in dihedral] - dihedral_matched = tuple([rename_dict[index] for index in dihedral]) - angle = rdkit_wrapper.get_dihedral(mol, *dihedral_matched, units="radians") - - pdb_dihedrals.append(names+[angle]) - self.dihedral_library[cluster_pdb] = pdb_dihedrals - - def save(self, output_path): - """ - Save the dihedral library - - Parameters - ---------- - output_path: str - Where to save the dihedral library - """ - with open(output_path, "w") as fw: - fw.write("* DIHEDRAL LIBRARY FILE\n") - fw.write('* File generated with peleffy-{}\n'.format( peleffy.__version__)) - for dihedral_file, dihedral_values in self.dihedral_library.items(): - fw.write("* File: {:s}\n".format(dihedral_file)) - fw.write("{:s} {:d} {:d}\n".format(self._molecule.tag.upper(), len(dihedral_values), len(self.dihedral_library))) - for dihedral in dihedral_values: - fw.write("{:s} {:s} {:s} {:s} {:5f}\n".format(*dihedral)) - fw.write("ENDDIHEDRALS\n") - fw.write("END\n") - - def list_all_dihedrals(self): - dihedrals = [] - seen_dihedrals = set() - for proper in self._topology.propers: - dihedral_indexes = (proper.atom1_idx, proper.atom2_idx, proper.atom3_idx, proper.atom4_idx) - if dihedral_indexes in seen_dihedrals: - continue - seen_dihedrals.add(dihedral_indexes) - dihedrals.append(list(dihedral_indexes)) - return dihedrals - - - def order_clusters_min_distances(self, clusters): - """ - Order a list of structures in a way that jumps to the next structure - involves the minimum change - - Parameters - ---------- - - clusters: list - Path to the structures - """ - rdkit_wrapper = RDKitToolkitWrapper() - distances = np.zeros((len(clusters), len(clusters))) - cluster_molecules = [] - for cluster in clusters: - # use the input molecule as template since the cluster structures - # probably will not have proper stereochemistry - cluster_molecules.append(molecule.Molecule(cluster, connectivity_template=self._molecule.rdkit_molecule)) - for i, cluster_mol in enumerate(cluster_molecules): - for j, cluster_mol_2 in enumerate(cluster_molecules[i+1:], start=i+1): - rmsd_value = rdkit_wrapper.get_rmsd(cluster_mol, cluster_mol_2) - distances[i, j] = rmsd_value - distances[j, i] = rmsd_value - return find_optimal_path_from_matrix(distances) - - -def find_optimal_path_from_matrix(distances): - """ - Find a Minimum Spanning Tree from a distance matrix. It uses networkx to - create the graph and the MST - - Parameters - ---------- - - Returns - ------- - min_path : list - Nodes on the graph that minimise heuristically the total distance - """ - graph = nx.from_numpy_matrix(distances) - min_dist = 1e6 - min_path = None - for node in graph: - dist, path = find_heuristic_path(graph.copy(), distances, node) - if dist < min_dist: - min_dist = dist - min_path = path - return min_path - -def find_heuristic_path(graph, distances, start_node): - n = distances.shape[0] - path = [start_node] - dist = 0 - while len(path) < n: - node = path[-1] - new_neighbor = min(graph[node].items(), key=lambda x: x[1]["weight"]) - graph.remove_node(node) - path.append(new_neighbor[0]) - dist += new_neighbor[1]['weight'] - dist += distances[path[-1], start_node] - return dist, path - -def read_dihedral_info_file(file_name): - """ - Read an ndx file containing the classification of all dihedrals and return - the flexible ones - - Parameters - ---------- - file_name: str - Path of the dihedrals file - - Returns - ------- - dihedrals : list - List of atom indices correponding to the flexible dihedrals - """ - read_dihedral = False - dihedrals = [] - with open(file_name) as f: - for line in f: - # reached the block of all dihedrals, we can stop reading - if line.startswith("[ ALL ]"): - return dihedrals - if line.startswith("["): - name = line.strip("\n[] ") - read_dihedral = name.startswith("F") - else: - if read_dihedral: - # the indices of the bce server are 1-based, while our - # indices are 0-based - dihedrals.append([int(x)-1 for x in line.rstrip().split()]) - read_dihedral = False - return dihedrals diff --git a/peleffy/BCEDihedrals/__init__.py b/peleffy/BCEDihedrals/__init__.py deleted file mode 100644 index 1d38f1af..00000000 --- a/peleffy/BCEDihedrals/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from .BCEDihedrals import BCEDihedrals diff --git a/peleffy/main.py b/peleffy/main.py index 95bc7ed8..6a5a28e6 100644 --- a/peleffy/main.py +++ b/peleffy/main.py @@ -55,13 +55,9 @@ def parse_args(args): parser.add_argument("-o", "--output", metavar="PATH", help="Output path. Default is the current working " + "directory") - parser.add_argument("--dihedrals_info_path", default=None, type=str, + parser.add_argument("--conformations_info_path", default=None, type=str, help="Path to the folder containing the BCE output" - + " used to collect dihedral angles for PELE") - parser.add_argument("--dihedrals_mode", default="all_dihedrals", type=str, - help="Mode for extracting the dihedrals: 'all_dihedrals' will"+ - "extract the values for all, while 'flexible_dihedrals' will extract"+ - "the values only for those that are marked as flexible") + + " used to collect conformations for PELE") parser.add_argument('--with_solvent', dest='with_solvent', help="Generate solvent parameters for OBC", action='store_true') @@ -107,7 +103,7 @@ def run_peleffy(pdb_file, charges_from_file=None, exclude_terminal_rotamers=True, output=None, with_solvent=False, as_datalocal=False, - dihedral_path=None, dihedral_mode="all_dihedrals"): + conformation_path=None): """ It runs peleffy. @@ -135,7 +131,7 @@ def run_peleffy(pdb_file, as_datalocal : bool Whether to save output files following PELE's DataLocal hierarchy or not - dihedral_path: str + conformation_path: str Path to the BCE server outupt to use to extract dihedral angles dihedral_mode: str Select what kind of dihedrals to extract (all or only flexible) @@ -167,7 +163,7 @@ def run_peleffy(pdb_file, from peleffy.topology import Molecule from peleffy.template import Impact from peleffy.solvent import OBC2 - from peleffy.BCEDihedrals import BCEDihedrals + from peleffy.BCEConformations import BCEConformations from peleffy.forcefield import ForceFieldSelector from peleffy.topology import Topology from peleffy.utils import parse_charges_from_mae @@ -187,8 +183,8 @@ def run_peleffy(pdb_file, output_path=output, as_datalocal=as_datalocal) - if dihedral_path is None: - # if dihedral_path is set, we don't want a rotamer library + if conformation_path is None: + # if conformation_path is set, we don't want a rotamer library rotamer_library = peleffy.topology.RotamerLibrary(molecule) rotamer_library.to_file(output_handler.get_rotamer_library_path()) @@ -219,19 +215,19 @@ def run_peleffy(pdb_file, solvent = OBC2(topology) solvent.to_file(output_handler.get_solvent_template_path()) - if dihedral_path is not None: + if conformation_path is not None: previous_level = log.get_level() - dihedrals = BCEDihedrals(topology, dihedral_path, mode=dihedral_mode) - dihedrals.calculate() - dihedrals.save(output_handler.get_dihedral_library_path()) + conformations = BCEConformations(topology, conformation_path) + conformations.calculate() + conformations.save(output_handler.get_conformation_library_path()) log.set_level(previous_level) log.info(' - All files were generated successfully:') - if dihedral_path is None: + if conformation_path is None: log.info(' - {}'.format(output_handler.get_rotamer_library_path())) log.info(' - {}'.format(output_handler.get_impact_template_path())) - if dihedral_path is not None: - log.info(' - {}'.format(output_handler.get_dihedral_library_path())) + if conformation_path is not None: + log.info(' - {}'.format(output_handler.get_conformation_library_path())) if with_solvent: log.info(' - {}'.format(output_handler.get_solvent_template_path())) @@ -280,9 +276,8 @@ def main(args): output=args.output, with_solvent=args.with_solvent, as_datalocal=args.as_datalocal, - dihedral_path=args.dihedrals_info_path, - charges_from_file=args.charges_from_file, - dihedral_mode=args.dihedrals_mode) + conformation_path=args.conformations_info_path, + charges_from_file=args.charges_from_file) if __name__ == '__main__': diff --git a/peleffy/utils/output.py b/peleffy/utils/output.py index 194eab65..8e206130 100644 --- a/peleffy/utils/output.py +++ b/peleffy/utils/output.py @@ -18,8 +18,8 @@ class OutputPathHandler(object): OPLS_IMPACT_TEMPLATE_PATH = 'DataLocal/Templates/OPLS2005/HeteroAtoms/' ROTAMER_LIBRARY_PATH = 'DataLocal/LigandRotamerLibs/' SOLVENT_TEMPLATE_PATH = 'DataLocal/OBC/' - DIHEDRALS_LIBRARY_PATH = 'DataLocal/Dihedrals/' - FILE_TYPES = ['impact template', 'rotamer library', 'solvent template', 'dihedral library'] + CONFORMATION_LIBRARY_PATH = 'DataLocal/Conformations/' + FILE_TYPES = ['impact template', 'rotamer library', 'solvent template', 'conformation library'] def __init__(self, molecule, forcefield, as_datalocal=False, output_path=None): @@ -67,7 +67,7 @@ def get_path(self, file_type, create_missing_folders=True): ---------- file_type : str The file type whose path is requested. One of - ['impact template', 'rotamer library', 'solvent template'] + ['impact template', 'rotamer library', 'solvent template', 'conformation library'] create_missing_folders : bool Whether to create missing folders or not. Default is True @@ -89,8 +89,8 @@ def get_path(self, file_type, create_missing_folders=True): if file_type.lower() == 'solvent template': return self.get_solvent_template_path(create_missing_folders) - if file_type.lower() == 'dihedral library': - return self.get_dihedral_library_path(create_missing_folders) + if file_type.lower() == 'conformation library': + return self.get_conformation_library_path(create_missing_folders) def get_impact_template_path(self, create_missing_folders=True): """ @@ -182,9 +182,9 @@ def get_solvent_template_path(self, create_missing_folders=True): return os.path.join(path, file_name) - def get_dihedral_library_path(self, create_missing_folders=True): + def get_conformation_library_path(self, create_missing_folders=True): """ - It returns the path for a dihedral library file. + It returns the path for a conformation library file. Parameters ---------- @@ -194,12 +194,12 @@ def get_dihedral_library_path(self, create_missing_folders=True): Returns ------- file_path : str - The path for a dihedral library file + The path for a conformation library file """ - file_name = self._molecule.tag.upper() + '.dihedral' + file_name = self._molecule.tag.upper() + '.conformation' if self.as_datalocal: - path = os.path.join(self.output_path, self.DIHEDRALS_LIBRARY_PATH) + path = os.path.join(self.output_path, self.CONFORMATION_LIBRARY_PATH) else: path = self.output_path From fc85440d046a8d64699f4f9e15a87e58aa6004b2 Mon Sep 17 00:00:00 2001 From: jgilaber Date: Mon, 15 Mar 2021 13:13:16 +0100 Subject: [PATCH 26/27] Update test for BCEConformations --- peleffy/data/parameters/ETH.conformation | 12 ++++ peleffy/tests/test_BCEConformations.py | 61 ++++++++++++++++++++ peleffy/tests/test_BCEDihedrals.py | 71 ------------------------ peleffy/tests/test_integrity.py | 2 +- peleffy/tests/test_main.py | 57 +++++++------------ peleffy/tests/test_utils.py | 64 ++++++++++----------- 6 files changed, 125 insertions(+), 142 deletions(-) create mode 100644 peleffy/data/parameters/ETH.conformation create mode 100644 peleffy/tests/test_BCEConformations.py delete mode 100644 peleffy/tests/test_BCEDihedrals.py diff --git a/peleffy/data/parameters/ETH.conformation b/peleffy/data/parameters/ETH.conformation new file mode 100644 index 00000000..1a7db46c --- /dev/null +++ b/peleffy/data/parameters/ETH.conformation @@ -0,0 +1,12 @@ +* CONFORMATION LIBRARY FILE +* File generated with peleffy-1.0.0+164.g06140f6 +* File: ../../PELEPipelineSanja/4S0V-LIG/CLUSTERS/CL1/cluster1.min.imaged.pdb +UNK 6 1 +_C1_ 0.000000 0.000000 0.000000 +_C2_ -1.305000 0.181000 -0.014000 +_H1_ 0.720000 0.822000 0.042000 +_H2_ 0.448000 -0.988000 -0.029000 +_H3_ -1.730000 1.171000 0.015000 +_H4_ -1.936000 -0.687000 -0.056000 +ENDCONFORMATION +END diff --git a/peleffy/tests/test_BCEConformations.py b/peleffy/tests/test_BCEConformations.py new file mode 100644 index 00000000..4a639091 --- /dev/null +++ b/peleffy/tests/test_BCEConformations.py @@ -0,0 +1,61 @@ +""" +This module contains tests that check that the BCEConformations module. +""" + +import numpy as np +import pytest + +def build_mock_BCEConformations(pdb_file): + from peleffy.topology import Molecule + from peleffy.BCEConformations import BCEConformations + from peleffy.forcefield import ForceFieldSelector + from peleffy.topology import Topology + + molecule = Molecule(pdb_file) + ff_selector = ForceFieldSelector() + forcefield = ff_selector.get_by_name("opls2005") + parameters = forcefield.parameterize(molecule, charge_method="opls2005") + topology = Topology(molecule, parameters) + return BCEConformations(topology, "") + +class TestBCEConformations(object): + """BCEConformations test.""" + + from peleffy.BCEConformations import BCEConformations + + def test_calculate_cluster_offsets(self): + from peleffy.utils import get_data_file_path + + pdb_path = get_data_file_path('ligands/ethylene.pdb') + golden_offsets = [["_C1_", 0.0, 0.0, 0.0], + ["_C2_", -1.305, 0.181, -0.014], + ["_H1_", 0.72, 0.822, 0.042], + ["_H2_", 0.448, -0.988, -0.029], + ["_H3_", -1.73, 1.171, 0.015], + ["_H4_", -1.936, -0.687, -0.056]] + bce_obj = build_mock_BCEConformations(pdb_path) + bce_obj.calculate_cluster_offsets(pdb_path) + for dih1, dih2 in zip(bce_obj.conformations_library[pdb_path], golden_offsets): + assert dih1[0] == dih2[0] + np.testing.assert_almost_equal(dih1[1:], dih2[1:], decimal=3) + + def test_write_dihedral_library(self): + import os + import tempfile + from peleffy.utils import get_data_file_path, temporary_cd + + pdb_path = get_data_file_path('ligands/ethylene.pdb') + bce_obj = build_mock_BCEConformations(pdb_path) + bce_obj.calculate_cluster_offsets(pdb_path) + calculated_lines = [] + with tempfile.TemporaryDirectory() as tmpdir: + with temporary_cd(tmpdir): + bce_obj.save(os.path.join(tmpdir, "ETH.conformation")) + with open(os.path.join(tmpdir, "ETH.conformation")) as f: + calculated_lines = f.readlines() + calculated_lines = calculated_lines[3:] + golden_conformation_library_path = get_data_file_path('parameters/ETH.conformation') + with open(golden_conformation_library_path) as f: + golden_lines = f.readlines()[3:] + + assert golden_lines == calculated_lines diff --git a/peleffy/tests/test_BCEDihedrals.py b/peleffy/tests/test_BCEDihedrals.py deleted file mode 100644 index c37052ed..00000000 --- a/peleffy/tests/test_BCEDihedrals.py +++ /dev/null @@ -1,71 +0,0 @@ -""" -This module contains tests that check that the BCEDihedrals module. -""" - -import numpy as np -import pytest - -def build_mock_BCEDihedrals(pdb_file): - from peleffy.topology import Molecule - from peleffy.BCEDihedrals import BCEDihedrals - from peleffy.forcefield import ForceFieldSelector - from peleffy.topology import Topology - - molecule = Molecule(pdb_file) - ff_selector = ForceFieldSelector() - forcefield = ff_selector.get_by_name("opls2005") - parameters = forcefield.parameterize(molecule, charge_method="opls2005") - topology = Topology(molecule, parameters) - return BCEDihedrals(topology, "", mode="") - -class TestBCEDihderals(object): - """BCEDihedrals test.""" - - from peleffy.BCEDihedrals import BCEDihedrals - - def test_list_dihedrals(self): - from peleffy.utils import get_data_file_path - - pdb_path = get_data_file_path('ligands/ethylene.pdb') - golden_dihedrals = [(2, 0, 1, 4), (2, 0, 1, 5), (3, 0, 1, 4), (3, 0, 1, 5)] - bce_obj = build_mock_BCEDihedrals(pdb_path) - dihedrals_indices = bce_obj.list_all_dihedrals() - dihedrals_indices.sort() - np.testing.assert_array_equal(dihedrals_indices, golden_dihedrals) - - def test_calculate_cluster_angles(self): - from peleffy.utils import get_data_file_path - - golden_dihedrals = [(2, 0, 1, 4), (2, 0, 1, 5), (3, 0, 1, 4), (3, 0, 1, 5)] - pdb_path = get_data_file_path('ligands/ethylene.pdb') - golden_angles = [ ["_H1_", "_C1_", "_C2_", "_H3_", 0.001151], - ["_H1_", "_C1_", "_C2_", "_H4_", -3.141278], - ["_H2_", "_C1_", "_C2_", "_H3_", -3.141278], - ["_H2_", "_C1_", "_C2_", "_H4_", -0.000366]] - bce_obj = build_mock_BCEDihedrals(pdb_path) - bce_obj.calculate_cluster_angles(pdb_path, golden_dihedrals, match_indexes=False) - for dih1, dih2 in zip(bce_obj.dihedral_library[pdb_path], golden_angles): - np.testing.assert_array_equal(dih1[:4], dih2[:4]) - np.testing.assert_almost_equal(dih1[-1], dih2[-1], decimal=3) - - def test_write_dihedral_library(self): - import os - import tempfile - from peleffy.utils import get_data_file_path, temporary_cd - - golden_dihedrals = [(2, 0, 1, 4), (2, 0, 1, 5), (3, 0, 1, 4), (3, 0, 1, 5)] - pdb_path = get_data_file_path('ligands/ethylene.pdb') - bce_obj = build_mock_BCEDihedrals(pdb_path) - bce_obj.calculate_cluster_angles(pdb_path, golden_dihedrals, match_indexes=False) - calculated_lines = [] - with tempfile.TemporaryDirectory() as tmpdir: - with temporary_cd(tmpdir): - bce_obj.save(os.path.join(tmpdir, "ETH.dihedral")) - with open(os.path.join(tmpdir, "ETH.dihedral")) as f: - calculated_lines = f.readlines() - calculated_lines = calculated_lines[3:] - golden_dihedral_library_path = get_data_file_path('parameters/ETH.dihedral') - with open(golden_dihedral_library_path) as f: - golden_lines = f.readlines()[3:] - - assert golden_lines == calculated_lines diff --git a/peleffy/tests/test_integrity.py b/peleffy/tests/test_integrity.py index f1871fd4..6b254649 100644 --- a/peleffy/tests/test_integrity.py +++ b/peleffy/tests/test_integrity.py @@ -17,7 +17,7 @@ def test_modules(self): from peleffy import template from peleffy import topology from peleffy import utils - from peleffy import BCEDihedrals + from peleffy import BCEConformations except ImportError as e: raise AssertionError('The following peleffy module is missing: ' diff --git a/peleffy/tests/test_main.py b/peleffy/tests/test_main.py index cd4b3ae8..3647c8d7 100644 --- a/peleffy/tests/test_main.py +++ b/peleffy/tests/test_main.py @@ -81,10 +81,8 @@ def test_peleffy_argparse(self): 'Unexpected silent settings were parsed' assert parsed_args.with_solvent is False, \ 'Unexpected with_solvent settings were parsed' - assert parsed_args.dihedrals_info_path is None, \ - 'Unexpected dihedral_path settings were passed' - assert parsed_args.dihedrals_mode == "all_dihedrals", \ - 'Unexpected dihedrals_mode settings were passed' + assert parsed_args.conformations_info_path is None, \ + 'Unexpected conformations_path settings were passed' # Test custom shorts parsed_args = parse_args(['toluene.pdb', @@ -113,10 +111,8 @@ def test_peleffy_argparse(self): 'Unexpected silent settings were parsed' assert parsed_args.with_solvent is False, \ 'Unexpected with_solvent settings were parsed' - assert parsed_args.dihedrals_info_path is None, \ - 'Unexpected dihedral_path settings were passed' - assert parsed_args.dihedrals_mode == "all_dihedrals", \ - 'Unexpected dihedrals_mode settings were passed' + assert parsed_args.conformations_info_path is None, \ + 'Unexpected conformations_path settings were passed' # Test custom longs parsed_args = parse_args(['methane.pdb', @@ -146,10 +142,8 @@ def test_peleffy_argparse(self): 'Unexpected silent settings were parsed' assert parsed_args.with_solvent is False, \ 'Unexpected with_solvent settings were parsed' - assert parsed_args.dihedrals_info_path is None, \ - 'Unexpected dihedral_path settings were passed' - assert parsed_args.dihedrals_mode == "all_dihedrals", \ - 'Unexpected dihedrals_mode settings were passed' + assert parsed_args.conformations_info_path is None, \ + 'Unexpected conformations_path settings were passed' # Test unexpected charge method with pytest.raises(SystemExit) as pytest_wrapped_e: @@ -182,10 +176,8 @@ def test_peleffy_argparse(self): 'Unexpected silent settings were parsed' assert parsed_args.with_solvent is False, \ 'Unexpected with_solvent settings were parsed' - assert parsed_args.dihedrals_info_path is None, \ - 'Unexpected dihedral_path settings were passed' - assert parsed_args.dihedrals_mode == "all_dihedrals", \ - 'Unexpected dihedrals_mode settings were passed' + assert parsed_args.conformations_info_path is None, \ + 'Unexpected conformations_path settings were passed' # Test charges_from_file argument parsed_args = parse_args(['BHP.pdb', @@ -213,10 +205,8 @@ def test_peleffy_argparse(self): 'Unexpected silent settings were parsed' assert parsed_args.with_solvent is False, \ 'Unexpected with_solvent settings were parsed' - assert parsed_args.dihedrals_info_path is None, \ - 'Unexpected dihedral_path settings were passed' - assert parsed_args.dihedrals_mode == "all_dihedrals", \ - 'Unexpected dihedrals_mode settings were passed' + assert parsed_args.conformations_info_path is None, \ + 'Unexpected conformations_path settings were passed' # Test include_terminal_rotamers argument parsed_args = parse_args(['methane.pdb', @@ -242,10 +232,8 @@ def test_peleffy_argparse(self): 'Unexpected silent settings were parsed' assert parsed_args.with_solvent is False, \ 'Unexpected with_solvent settings were parsed' - assert parsed_args.dihedrals_info_path is None, \ - 'Unexpected dihedral_path settings were passed' - assert parsed_args.dihedrals_mode == "all_dihedrals", \ - 'Unexpected dihedrals_mode settings were passed' + assert parsed_args.conformations_info_path is None, \ + 'Unexpected conformations_path settings were passed' # Test silent argument parsed_args = parse_args(['methane.pdb', @@ -271,10 +259,8 @@ def test_peleffy_argparse(self): 'Unexpected silent settings were parsed' assert parsed_args.with_solvent is False, \ 'Unexpected with_solvent settings were parsed' - assert parsed_args.dihedrals_info_path is None, \ - 'Unexpected dihedral_path settings were passed' - assert parsed_args.dihedrals_mode == "all_dihedrals", \ - 'Unexpected dihedrals_mode settings were passed' + assert parsed_args.conformations_info_path is None, \ + 'Unexpected conformations_path settings were passed' parse_args(['methane.pdb', '-s']) == parse_args(['methane.pdb', '--silent']) @@ -303,18 +289,15 @@ def test_peleffy_argparse(self): 'Unexpected silent settings were parsed' assert parsed_args.with_solvent is False, \ 'Unexpected with_solvent settings were parsed' - assert parsed_args.dihedrals_info_path is None, \ - 'Unexpected dihedral_path settings were passed' - assert parsed_args.dihedrals_mode == "all_dihedrals", \ - 'Unexpected dihedrals_mode settings were passed' + assert parsed_args.conformations_info_path is None, \ + 'Unexpected conformations_path settings were passed' parse_args(['methane.pdb', '-d']) == parse_args(['methane.pdb', '--debug']) # Test dihedral library arguments parsed_args = parse_args(['methane.pdb', - '--dihedrals_info_path', 'test_path', - '--dihedrals_mode', 'flexible_dihedrals']) + '--conformations_info_path', 'test_path']) assert parsed_args.as_datalocal is False, \ 'Unexpected as_datalocal settings were parsed' @@ -336,10 +319,8 @@ def test_peleffy_argparse(self): 'Unexpected silent settings were parsed' assert parsed_args.with_solvent is False, \ 'Unexpected with_solvent settings were parsed' - assert parsed_args.dihedrals_info_path == "test_path", \ - 'Unexpected dihedral_path settings were passed' - assert parsed_args.dihedrals_mode == "flexible_dihedrals", \ - 'Unexpected dihedrals_mode settings were passed' + assert parsed_args.conformations_info_path == "test_path", \ + 'Unexpected conformations_path settings were passed' def test_peleffy_main(self): """It checks the main function of peleffy.""" diff --git a/peleffy/tests/test_utils.py b/peleffy/tests/test_utils.py index 028e69f4..62f9f268 100644 --- a/peleffy/tests/test_utils.py +++ b/peleffy/tests/test_utils.py @@ -194,9 +194,9 @@ def test_non_datalocal_paths(self): assert output_handler.get_solvent_template_path() == \ './ligandParams.txt', \ 'Unexpected default solvent parameters path' - assert output_handler.get_dihedral_library_path() == \ - './{}.dihedral'.format(tag.upper()), \ - 'Unexpected default dihedral library path' + assert output_handler.get_conformation_library_path() == \ + './{}.conformation'.format(tag.upper()), \ + 'Unexpected default conformation library path' # Initialize output handler with an output_path set with tempfile.TemporaryDirectory() as tmpdir: @@ -214,9 +214,9 @@ def test_non_datalocal_paths(self): assert output_handler.get_solvent_template_path() == \ os.path.join(tmpdir, 'output', 'ligandParams.txt'), \ 'Unexpected default solvent parameters path' - assert output_handler.get_dihedral_library_path() == \ - os.path.join(tmpdir, 'output', '{}.dihedral'.format(tag.upper())), \ - 'Unexpected default dihedral library path' + assert output_handler.get_conformation_library_path() == \ + os.path.join(tmpdir, 'output', '{}.conformation'.format(tag.upper())), \ + 'Unexpected default conformation library path' def test_datalocal_paths_for_openff(self): """It tests the datalocal paths assignment for OpenFF.""" @@ -250,10 +250,10 @@ def test_datalocal_paths_for_openff(self): create_missing_folders=False) == \ './DataLocal/OBC/ligandParams.txt', \ 'Unexpected default solvent parameters path' - assert output_handler.get_dihedral_library_path( + assert output_handler.get_conformation_library_path( create_missing_folders=False) == \ - './DataLocal/Dihedrals/'+'{}.dihedral'.format(tag.upper()), \ - 'Unexpected default dihedral library path' + './DataLocal/Conformations/'+'{}.conformation'.format(tag.upper()), \ + 'Unexpected default conformation library path' with tempfile.TemporaryDirectory() as tmpdir: output_handler = OutputPathHandler( @@ -275,10 +275,10 @@ def test_datalocal_paths_for_openff(self): create_missing_folders=False) == \ tmpdir + '/output/DataLocal/OBC/ligandParams.txt', \ 'Unexpected default solvent parameters path' - assert output_handler.get_dihedral_library_path( + assert output_handler.get_conformation_library_path( create_missing_folders=False) == \ - tmpdir + '/output/DataLocal/Dihedrals/'+'{}.dihedral'.format(tag.upper()), \ - 'Unexpected default dihedral library path' + tmpdir + '/output/DataLocal/Conformations/'+'{}.conformation'.format(tag.upper()), \ + 'Unexpected default conformation library path' def test_datalocal_paths_for_opls(self): """It tests the datalocal paths assignment for OPLS2005.""" @@ -313,10 +313,10 @@ def test_datalocal_paths_for_opls(self): create_missing_folders=False) == \ './DataLocal/OBC/ligandParams.txt', \ 'Unexpected default solvent parameters path' - assert output_handler.get_dihedral_library_path( + assert output_handler.get_conformation_library_path( create_missing_folders=False) == \ - './DataLocal/Dihedrals/'+'{}.dihedral'.format(tag.upper()), \ - 'Unexpected default dihedral library path' + './DataLocal/Conformations/'+'{}.conformation'.format(tag.upper()), \ + 'Unexpected default conformation library path' # Initialize output handler with an output_path set with tempfile.TemporaryDirectory() as tmpdir: @@ -339,10 +339,10 @@ def test_datalocal_paths_for_opls(self): os.path.join(tmpdir, 'output', 'DataLocal/OBC/ligandParams.txt'), \ 'Unexpected default solvent parameters path' - assert output_handler.get_dihedral_library_path( + assert output_handler.get_conformation_library_path( create_missing_folders=False) == \ - os.path.join(tmpdir, 'output/DataLocal/Dihedrals/'+'{}.dihedral'.format(tag.upper())), \ - 'Unexpected default dihedral library path' + os.path.join(tmpdir, 'output/DataLocal/Conformations/'+'{}.conformation'.format(tag.upper())), \ + 'Unexpected default conformation library path' def test_datalocal_paths_for_offopls(self): """ @@ -382,10 +382,10 @@ def test_datalocal_paths_for_offopls(self): create_missing_folders=False) == \ './DataLocal/OBC/ligandParams.txt', \ 'Unexpected default solvent parameters path' - assert output_handler.get_dihedral_library_path( + assert output_handler.get_conformation_library_path( create_missing_folders=False) == \ - './DataLocal/Dihedrals/'+'{}.dihedral'.format(tag.upper()), \ - 'Unexpected default dihedral library path' + './DataLocal/Conformations/'+'{}.conformation'.format(tag.upper()), \ + 'Unexpected default conformation library path' # Initialize output handler with an output_path set with tempfile.TemporaryDirectory() as tmpdir: @@ -408,10 +408,10 @@ def test_datalocal_paths_for_offopls(self): os.path.join(tmpdir, 'output', 'DataLocal/OBC/ligandParams.txt'), \ 'Unexpected default solvent parameters path' - assert output_handler.get_dihedral_library_path( + assert output_handler.get_conformation_library_path( create_missing_folders=False) == \ - os.path.join(tmpdir, 'output/DataLocal/Dihedrals/'+'{}.dihedral'.format(tag.upper())), \ - 'Unexpected default dihedral library path' + os.path.join(tmpdir, 'output/DataLocal/Conformations/'+'{}.conformation'.format(tag.upper())), \ + 'Unexpected default conformation library path' # Initialize output handler without output_path output_handler = OutputPathHandler(molecule, hybridff, @@ -435,10 +435,10 @@ def test_datalocal_paths_for_offopls(self): create_missing_folders=False) == \ './DataLocal/OBC/ligandParams.txt', \ 'Unexpected default solvent parameters path' - assert output_handler.get_dihedral_library_path( + assert output_handler.get_conformation_library_path( create_missing_folders=False) == \ - './DataLocal/Dihedrals/'+'{}.dihedral'.format(tag.upper()), \ - 'Unexpected default dihedral library path' + './DataLocal/Conformations/'+'{}.conformation'.format(tag.upper()), \ + 'Unexpected default conformation library path' # Initialize output handler with an output_path set with tempfile.TemporaryDirectory() as tmpdir: @@ -461,11 +461,11 @@ def test_datalocal_paths_for_offopls(self): os.path.join(tmpdir, 'output', 'DataLocal/OBC/ligandParams.txt'), \ 'Unexpected default solvent parameters path' - assert output_handler.get_dihedral_library_path( + assert output_handler.get_conformation_library_path( create_missing_folders=False) == \ os.path.join(tmpdir, 'output', - 'DataLocal/Dihedrals/'+'{}.dihedral'.format(tag.upper())), \ - 'Unexpected default dihedral library path' + 'DataLocal/Conformations/'+'{}.conformation'.format(tag.upper())), \ + 'Unexpected default conformation library path' def test_folder_creation(self): """ @@ -503,7 +503,7 @@ def test_folder_creation(self): path_dir = os.path.dirname(path) assert os.path.isdir(path_dir) is False, \ 'This directory should not exist' - path = output_handler.get_dihedral_library_path( + path = output_handler.get_conformation_library_path( create_missing_folders=False) path_dir = os.path.dirname(path) assert os.path.isdir(path_dir) is False, \ @@ -533,7 +533,7 @@ def test_folder_creation(self): path_dir = os.path.dirname(path) assert os.path.isdir(path_dir) is True, \ 'This directory should exist' - path = output_handler.get_dihedral_library_path( + path = output_handler.get_conformation_library_path( create_missing_folders=True) path_dir = os.path.dirname(path) assert os.path.isdir(path_dir) is True, \ From a074110fccf4b315f7968eed3a8fea405e0a60e2 Mon Sep 17 00:00:00 2001 From: jgilaber Date: Wed, 7 Apr 2021 12:24:58 +0200 Subject: [PATCH 27/27] Use parameterize_opls2005 function to avoid errors --- peleffy/tests/test_BCEConformations.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/peleffy/tests/test_BCEConformations.py b/peleffy/tests/test_BCEConformations.py index 4a639091..4b9aa65c 100644 --- a/peleffy/tests/test_BCEConformations.py +++ b/peleffy/tests/test_BCEConformations.py @@ -5,16 +5,17 @@ import numpy as np import pytest -def build_mock_BCEConformations(pdb_file): +def build_mock_BCEConformations(pdb_file, ffld_file): from peleffy.topology import Molecule from peleffy.BCEConformations import BCEConformations from peleffy.forcefield import ForceFieldSelector from peleffy.topology import Topology + from .utils import parameterize_opls2005 molecule = Molecule(pdb_file) ff_selector = ForceFieldSelector() forcefield = ff_selector.get_by_name("opls2005") - parameters = forcefield.parameterize(molecule, charge_method="opls2005") + parameters = parameterize_opls2005(forcefield, molecule, ffld_file, charge_method="opls2005") topology = Topology(molecule, parameters) return BCEConformations(topology, "") @@ -27,13 +28,14 @@ def test_calculate_cluster_offsets(self): from peleffy.utils import get_data_file_path pdb_path = get_data_file_path('ligands/ethylene.pdb') + ffld_path = get_data_file_path("tests/ETL_ffld_output.txt") golden_offsets = [["_C1_", 0.0, 0.0, 0.0], ["_C2_", -1.305, 0.181, -0.014], ["_H1_", 0.72, 0.822, 0.042], ["_H2_", 0.448, -0.988, -0.029], ["_H3_", -1.73, 1.171, 0.015], ["_H4_", -1.936, -0.687, -0.056]] - bce_obj = build_mock_BCEConformations(pdb_path) + bce_obj = build_mock_BCEConformations(pdb_path, ffld_path) bce_obj.calculate_cluster_offsets(pdb_path) for dih1, dih2 in zip(bce_obj.conformations_library[pdb_path], golden_offsets): assert dih1[0] == dih2[0] @@ -45,7 +47,8 @@ def test_write_dihedral_library(self): from peleffy.utils import get_data_file_path, temporary_cd pdb_path = get_data_file_path('ligands/ethylene.pdb') - bce_obj = build_mock_BCEConformations(pdb_path) + ffld_path = get_data_file_path("tests/ETL_ffld_output.txt") + bce_obj = build_mock_BCEConformations(pdb_path, ffld_path) bce_obj.calculate_cluster_offsets(pdb_path) calculated_lines = [] with tempfile.TemporaryDirectory() as tmpdir: