diff --git a/pdb2pqr/pdb2pqr/ligand/mol2.py b/pdb2pqr/pdb2pqr/ligand/mol2.py index 9b5cbf1a2..533e275c1 100644 --- a/pdb2pqr/pdb2pqr/ligand/mol2.py +++ b/pdb2pqr/pdb2pqr/ligand/mol2.py @@ -5,6 +5,7 @@ """ import logging from collections import OrderedDict +from itertools import combinations import numpy @@ -137,15 +138,79 @@ class Mol2Molecule: for torsion in atom.torsions: self.torsions.add(torsion) + @staticmethod + def rotate_to_smallest(path): + """Rotate cycle path so that it begins with the smallest node. + + This was borrowed from StackOverflow: https://j.mp/2AHaukj + + Args: + path: list of atom names + Returns: + rotated path (list) + """ + n = path.index(min(path)) + return path[n:]+path[:n] + + def find_new_rings(self, path, rings, level=0): + """Find new rings in molecule. + + This was borrowed from StackOverflow: https://j.mp/2AHaukj + + Args: + path: list of atom names + rings: current list of rings + level: recursion level + Returns: + new list of rings + """ + start_node = path[0] + next_node = None + sub_path = [] + for bond in self.bonds: + atom1 = bond.atoms[0] + atom2 = bond.atoms[1] + if start_node in (atom1, atom2): + if atom1 == start_node: + next_node = atom2 + else: + next_node = atom1 + if next_node not in path: + sub_path = [next_node] + sub_path.extend(path) + rings = self.find_new_rings(sub_path, rings, level+1) + elif len(path) > 2 and next_node == path[-1]: + path_ = self.rotate_to_smallest(path) + inv_path = tuple(self.rotate_to_smallest(path_[::-1])) + path_ = tuple(path_) + if (path_ not in rings) and (inv_path not in rings): + rings.add(tuple(path_)) + return rings + def set_rings(self): """Set all rings in molecule. - Like many things, this was borrowed from StackOverflow: - https://stackoverflow.com/questions/12367801/finding-all-cycles-in-undirected-graphs + This was borrowed from StackOverflow: https://j.mp/2AHaukj """ + self.rings = set() + rings = set() + # Generate all rings for bond in self.bonds: - _LOGGER.error(str(bond)) - raise NotImplementedError() + for atom_name in bond.atoms: + rings = self.find_new_rings([atom_name], rings) + # Prune rings that are products of other rings + ring_sets = [] + for i in range(2, len(rings)+1): + for combo in combinations(rings, i): + ring_set = set().union(*combo) + ring_sets.append(ring_set) + for ring in rings: + ring_set = set(ring) + if ring_set in ring_sets: + _LOGGER.debug("Fused ring: %s", ring) + else: + _LOGGER.debug("Unfused ring: %s", ring) + self.rings.add(ring) def read(self, mol2_file): """Routines for reading MOL2 file. @@ -254,5 +319,5 @@ class Mol2Molecule: atom2.bonded_atoms.append(atom1) self.bonds.append(bond) self.set_torsions() - # self.set_rings() + self.set_rings() return mol2_file diff --git a/pdb2pqr/tests/ligand_test.py b/pdb2pqr/tests/ligand_test.py index 89e3069a7..4e6f3820f 100644 --- a/pdb2pqr/tests/ligand_test.py +++ b/pdb2pqr/tests/ligand_test.py @@ -58,17 +58,53 @@ def test_torsions(input_mol2): mol2_path = Path("tests/data") / input_mol2 with open(mol2_path, "rt") as mol2_file: ligand.read(mol2_file) - try: - benchmark = TORSION_RESULTS[input_mol2] - diff = ligand.torsions ^ benchmark - if len(diff) > 0: - err = "Torsion test failed for %s: %s" % ( - input_mol2, sorted(list(diff))) - raise ValueError(err) - except KeyError: - _LOGGER.warning( - "Skipping torsions for %s: %s", input_mol2, - sorted(list(ligand.torsions))) + try: + benchmark = TORSION_RESULTS[input_mol2] + diff = ligand.torsions ^ benchmark + if len(diff) > 0: + err = "Torsion test failed for %s: %s" % ( + input_mol2, sorted(list(diff))) + raise ValueError(err) + except KeyError: + _LOGGER.warning( + "Skipping torsion test for %s: %s", input_mol2, + sorted(list(ligand.torsions))) + + +RING_RESULTS = { + "ethanol.mol2": set(), + "glycerol.mol2": set(), + "cyclohexane.mol2": {('CAA', 'CAD', 'CAE', 'CAF', 'CAC', 'CAB')}, + "naphthalene.mol2": { + ('CAA', 'CAB', 'CAC', 'CAH', 'CAG', 'CAF'), + ('CAC', 'CAH', 'CAI', 'CAJ', 'CAE', 'CAD')}, + "anthracene.mol2": { + ('CAC', 'CAJ', 'CAK', 'CAL', 'CAE', 'CAD'), + ('CAE', 'CAL', 'CAM', 'CAN', 'CAG', 'CAF'), + ('CAA', 'CAB', 'CAC', 'CAJ', 'CAI', 'CAH') + } +} + +@pytest.mark.parametrize("input_mol2", [ + "cyclohexane.mol2", "ethanol.mol2", "glycerol.mol2", "anthracene.mol2", + "naphthalene.mol2"]) +def test_rings(input_mol2): + """Test assignment of torsion angles.""" + ligand = parameterize.ParameterizedMolecule() + mol2_path = Path("tests/data") / input_mol2 + with open(mol2_path, "rt") as mol2_file: + ligand.read(mol2_file) + try: + benchmark = RING_RESULTS[input_mol2] + diff = ligand.rings ^ benchmark + if len(diff) > 0: + err = "Ring test failed for %s: %s" % ( + input_mol2, sorted(list(diff))) + raise ValueError(err) + except KeyError: + _LOGGER.warning( + "Skipping ring test for %s: %s", input_mol2, + sorted(list(ligand.rings))) @pytest.mark.parametrize("input_pdb", ["1HPX", "1QBS", "1US0"], ids=str)