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/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/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/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/main.py b/peleffy/main.py index 909be3a7..6a5a28e6 100644 --- a/peleffy/main.py +++ b/peleffy/main.py @@ -55,6 +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("--conformations_info_path", default=None, type=str, + help="Path to the folder containing the BCE output" + + " used to collect conformations for PELE") parser.add_argument('--with_solvent', dest='with_solvent', help="Generate solvent parameters for OBC", action='store_true') @@ -99,8 +102,8 @@ def run_peleffy(pdb_file, charge_method=DEFAULT_CHARGE_METHOD, charges_from_file=None, exclude_terminal_rotamers=True, - output=None, with_solvent=False, - as_datalocal=False): + output=None, with_solvent=False, as_datalocal=False, + conformation_path=None): """ It runs peleffy. @@ -128,6 +131,10 @@ def run_peleffy(pdb_file, as_datalocal : bool Whether to save output files following PELE's DataLocal hierarchy or not + 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) """ if charges_from_file is not None: charge_method_str = 'file\n' \ @@ -156,6 +163,7 @@ def run_peleffy(pdb_file, from peleffy.topology import Molecule from peleffy.template import Impact from peleffy.solvent import OBC2 + from peleffy.BCEConformations import BCEConformations from peleffy.forcefield import ForceFieldSelector from peleffy.topology import Topology from peleffy.utils import parse_charges_from_mae @@ -175,8 +183,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 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()) # Parameterize molecule with the selected force field log.info(' - Parameterizing molecule') @@ -205,9 +215,19 @@ def run_peleffy(pdb_file, solvent = OBC2(topology) solvent.to_file(output_handler.get_solvent_template_path()) + if conformation_path is not None: + previous_level = log.get_level() + 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:') - log.info(' - {}'.format(output_handler.get_rotamer_library_path())) + if conformation_path is None: + log.info(' - {}'.format(output_handler.get_rotamer_library_path())) log.info(' - {}'.format(output_handler.get_impact_template_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())) @@ -256,6 +276,7 @@ def main(args): output=args.output, with_solvent=args.with_solvent, as_datalocal=args.as_datalocal, + conformation_path=args.conformations_info_path, charges_from_file=args.charges_from_file) diff --git a/peleffy/tests/test_BCEConformations.py b/peleffy/tests/test_BCEConformations.py new file mode 100644 index 00000000..4b9aa65c --- /dev/null +++ b/peleffy/tests/test_BCEConformations.py @@ -0,0 +1,64 @@ +""" +This module contains tests that check that the BCEConformations module. +""" + +import numpy as np +import pytest + +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 = parameterize_opls2005(forcefield, molecule, ffld_file, 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') + 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, 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] + 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') + 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: + 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_integrity.py b/peleffy/tests/test_integrity.py index 475dad89..6b254649 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 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 b6e2e57f..3647c8d7 100644 --- a/peleffy/tests/test_main.py +++ b/peleffy/tests/test_main.py @@ -81,6 +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.conformations_info_path is None, \ + 'Unexpected conformations_path settings were passed' # Test custom shorts parsed_args = parse_args(['toluene.pdb', @@ -109,6 +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.conformations_info_path is None, \ + 'Unexpected conformations_path settings were passed' # Test custom longs parsed_args = parse_args(['methane.pdb', @@ -138,6 +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.conformations_info_path is None, \ + 'Unexpected conformations_path settings were passed' # Test unexpected charge method with pytest.raises(SystemExit) as pytest_wrapped_e: @@ -170,6 +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.conformations_info_path is None, \ + 'Unexpected conformations_path settings were passed' # Test charges_from_file argument parsed_args = parse_args(['BHP.pdb', @@ -197,6 +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.conformations_info_path is None, \ + 'Unexpected conformations_path settings were passed' # Test include_terminal_rotamers argument parsed_args = parse_args(['methane.pdb', @@ -222,6 +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.conformations_info_path is None, \ + 'Unexpected conformations_path settings were passed' # Test silent argument parsed_args = parse_args(['methane.pdb', @@ -247,6 +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.conformations_info_path is None, \ + 'Unexpected conformations_path settings were passed' parse_args(['methane.pdb', '-s']) == parse_args(['methane.pdb', '--silent']) @@ -275,10 +289,39 @@ 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.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', + '--conformations_info_path', 'test_path']) + + 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.conformations_info_path == "test_path", \ + 'Unexpected conformations_path settings were passed' + def test_peleffy_main(self): """It checks the main function of peleffy.""" from peleffy.main import parse_args, main 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) diff --git a/peleffy/tests/test_utils.py b/peleffy/tests/test_utils.py index d0705a59..62f9f268 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_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: @@ -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_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.""" @@ -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_conformation_library_path( + create_missing_folders=False) == \ + './DataLocal/Conformations/'+'{}.conformation'.format(tag.upper()), \ + 'Unexpected default conformation 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_conformation_library_path( + create_missing_folders=False) == \ + 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.""" @@ -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_conformation_library_path( + create_missing_folders=False) == \ + './DataLocal/Conformations/'+'{}.conformation'.format(tag.upper()), \ + 'Unexpected default conformation 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_conformation_library_path( + create_missing_folders=False) == \ + os.path.join(tmpdir, 'output/DataLocal/Conformations/'+'{}.conformation'.format(tag.upper())), \ + 'Unexpected default conformation 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_conformation_library_path( + create_missing_folders=False) == \ + './DataLocal/Conformations/'+'{}.conformation'.format(tag.upper()), \ + 'Unexpected default conformation 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_conformation_library_path( + create_missing_folders=False) == \ + 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, @@ -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_conformation_library_path( + create_missing_folders=False) == \ + './DataLocal/Conformations/'+'{}.conformation'.format(tag.upper()), \ + 'Unexpected default conformation 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_conformation_library_path( + create_missing_folders=False) == \ + os.path.join(tmpdir, 'output', + 'DataLocal/Conformations/'+'{}.conformation'.format(tag.upper())), \ + 'Unexpected default conformation 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_conformation_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_conformation_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): """ diff --git a/peleffy/utils/output.py b/peleffy/utils/output.py index 2359f158..8e206130 100644 --- a/peleffy/utils/output.py +++ b/peleffy/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'] + 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): @@ -66,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 @@ -88,6 +89,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() == 'conformation library': + return self.get_conformation_library_path(create_missing_folders) + def get_impact_template_path(self, create_missing_folders=True): """ It returns the path for an Impact template file. @@ -178,6 +182,32 @@ def get_solvent_template_path(self, create_missing_folders=True): return os.path.join(path, file_name) + def get_conformation_library_path(self, create_missing_folders=True): + """ + It returns the path for a conformation 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 conformation library file + """ + file_name = self._molecule.tag.upper() + '.conformation' + + if self.as_datalocal: + path = os.path.join(self.output_path, self.CONFORMATION_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): """ diff --git a/peleffy/utils/toolkits.py b/peleffy/utils/toolkits.py index 3b80684c..fd82efd9 100644 --- a/peleffy/utils/toolkits.py +++ b/peleffy/utils/toolkits.py @@ -282,6 +282,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 get_atom_degrees(self, molecule): """ It returns the ordered list of atom degrees. The degree of an atom @@ -527,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()): diff --git a/peleffy/utils/utils.py b/peleffy/utils/utils.py index 603a3f95..1edf0374 100644 --- a/peleffy/utils/utils.py +++ b/peleffy/utils/utils.py @@ -351,6 +351,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.