Use a list of oemols for PointMutationExecutor input instead of files #746

glass-w · 2020-10-26T11:28:52Z

This PR adapts the input of PointMutationExecutor to allow for a list of oemol objects to be used instead. Both @schallerdavid and I were having difficulty parsing ligand .sdf files containing a nitrogen that was perceived as needing to be tetrahedral at this line in perses.utils.openeye. This PR allows the user to bypass internal oemol preparation.

Changed:

PointMutationExecutor now takes inputs: .sdf, .pdb, [oemol].
ligand_file renamed to ligand_input to be more general.

An example script of how to prepare the [oemol] object is below, input names changed from my use case:

from perses.app.relative_point_mutation_setup import PointMutationExecutor
from perses.utils.openeye import generate_unique_atom_names
import oechem
from simtk import unit

ifs = oechem.oemolistream()
ifs.open('./ligand.sdf')
mol_list = [oechem.OEMol(mol) for mol in ifs.GetOEMols()]
mol = mol_list[0]

if len([atom.GetName() for atom in mol.GetAtoms()]) > len(set([atom.GetName() for atom in mol.GetAtoms()])):
        mol = generate_unique_atom_names(mol)
oechem.OEAssignAromaticFlags(mol, oechem.OEAroModelOpenEye)
oechem.OEAssignHybridization(mol)
oechem.OEAddExplicitHydrogens(mol)
oechem.OEPerceiveChiral(mol)

mol_list = [mol]
protein_file = './protein.pdb'

x = PointMutationExecutor(
protein_file,
'1',
'510', # an example residue number
'VAL', # an example mutation
ligand_input=mol_list,
ionic_strength=0.15 * unit.molar,
small_molecule_forcefields='openff-1.2.0'
)

codecov · 2020-10-26T11:48:11Z

Codecov Report

Merging #746 into master will increase coverage by 0.02%.
The diff coverage is n/a.

zhang-ivy

@dominicrufa can you check my comments below to make sure they are correct? after that, @glass-w can make the changes

zhang-ivy · 2020-10-26T19:02:51Z

perses/app/relative_point_mutation_setup.py

+                    ligand_n_atoms = ligand_md_topology.n_atoms
+
+            elif isinstance(ligand_input, list): # list of oemol object(s)
+                ligand_mol = ligand_input[0] # assume only one molecule in the list


If we're assuming there is only molecule in the list, let's just allow the options for ligand_input to be a .pdb, .sdf or oemol object (passing a list thats always one elements feels unnecessary). You should be able to check if the ligand_input is an instance of oechem.OEMol

don't know if reading a small molecule as a .pdb is safe, come to think of it

the small molecule should not be passed as a pdb. if ligand_input is a protein, it should be passed as a pdb. this will be noted in function description (it's in one of my suggestions in this PR)

zhang-ivy · 2020-10-26T19:04:05Z

perses/app/relative_point_mutation_setup.py

+                    ligand_n_atoms = ligand_md_topology.n_atoms
+
+            elif isinstance(ligand_input, list): # list of oemol object(s)
+                ligand_mol = ligand_input[0] # assume only one molecule in the list
                molecules.append(Molecule.from_openeye(ligand_mol, allow_undefined_stereo=False))


I think ligand_mol = generate_unique_atom_names(ligand_mol) is necessary to include before this line.. did you see any errors when you ran with this version of PointMutationExecutor?

oh wait, i just noticed generate_unique_atom_names is redundant.
it's called here:

perses/perses/utils/openeye.py

Line 354 in f640bf8

molecule = generate_unique_atom_names(molecule)

so it doesn't need to be called again after createMolFromSDF

zhang-ivy · 2020-10-26T19:06:04Z

perses/app/relative_point_mutation_setup.py

@@ -116,7 +116,7 @@ def __init__(self,
                if phase == vacuum, then the complex will not be solvated with water; else, it will be solvated with tip3p
            conduct_endstate_validation : bool, default True
                whether to conduct an endstate validation of the hybrid topology factory
-            ligand_file : str, default None
+            ligand_input : str or list of oemol objects, default None
                path to ligand of interest (i.e. small molecule or protein); .sdf or .pdb


can you change this to:

path to ligand of interest: .pdb for protein and .sdf or oemol for small molecule

zhang-ivy · 2020-10-26T19:11:34Z

perses/app/relative_point_mutation_setup.py

+            if isinstance(ligand_input, str):
+                if ligand_input.endswith('.sdf'): # small molecule
+                        ligand_mol = createOEMolFromSDF(ligand_input, index=ligand_index)
+                        ligand_mol = generate_unique_atom_names(ligand_mol)


This line is redundant, as `generate_unique_atom_names is called here:

perses/perses/utils/openeye.py

Line 354 in f640bf8

molecule = generate_unique_atom_names(molecule)

@dominicrufa can we remove it?

yep. remove in the relative_point_mutation_setup.py, you mean?

yes in relative_point_mutation_setup.py

zhang-ivy · 2020-10-26T19:18:35Z

@glass-w I noticed your script above that generates the oemol does not include the functionality covered by these lines in createOEMolFromSDF() ... Did you get any errors when running the above script? if not, it may be fine to not include the lines I linked (as long as when you inspect the old and new pdbs from the htf, the small molecule geometries look ok). but it is probably best practice to incorporate the lines i linked above in your script for generating the oemol (you'll just have to find a workaround that checks that carbons have tetrahedral stereochemistry and nitrogens are pyramidal)

jchodera · 2020-10-26T22:36:54Z

This PR adapts the input of PointMutationExecutor to allow for a list of oemol objects to be used instead.

Instead of OEMols, can we use openforcefield.topology.Molecule objects?
These will be much more portable in the future.

Internally, we can use molecule.to_oemol() whenever we need an OEMol().

dominicrufa · 2020-10-26T22:51:16Z

This PR adapts the input of PointMutationExecutor to allow for a list of oemol objects to be used instead.

Instead of OEMols, can we use openforcefield.topology.Molecule objects?
These will be much more portable in the future.

Internally, we can use molecule.to_oemol() whenever we need an OEMol().

is there anything lossy between oemols and OFF molecule objects? if not, then this is perfectly fine and can be added here.

there probably should be a test written somewhere here.

jchodera · 2020-10-27T01:40:25Z

is there anything lossy between oemols and OFF molecule objects?

I believe some non-chemical information, such as SD tags, would be lost. But atom ordering and geometry should be preserved.

Is there any particular information you're worried about being lost, @dominicrufa? If so, we can check (or you can just ask on the OFF slack).

glass-w · 2020-10-27T12:12:57Z

@zhang-ivy

Did you get any errors when running the above script?

I didn't get any errors when running the above script and the PointMutationExecutor successfully generated a series of htfs.

as long as when you inspect the old and new pdbs from the htf, the small molecule geometries look ok

Yes the geometries looked fine when visualising.

you'll just have to find a workaround that checks that carbons have tetrahedral stereochemistry and nitrogens are pyramidal

Yes this is a key issue when using OE3DToInternalStereo. It seems that for a subset of these ligands running OE3DToInternalStereo considers a nitrogen to be a stereocenter so then running has_undefined_stereocenters returns True. At the moment there doesn't seem to be a workaround (cc @schallerdavid).

I agree that if a user were to prepare their own OEMol object for parsing to PointMutationExecutor then they should apply both checks.

glass-w · 2020-10-27T12:20:21Z

@jchodera

Instead of OEMols, can we use openforcefield.topology.Molecule objects?

Yes I think this would be an idea for the future. Perhaps we can raise this as an issue and generate a separate PR? Right now this change is needed for functionality with this set of ligands.

zhang-ivy · 2020-10-27T13:45:47Z

perses/app/relative_point_mutation_setup.py

@@ -116,8 +116,8 @@ def __init__(self,
                if phase == vacuum, then the complex will not be solvated with water; else, it will be solvated with tip3p
            conduct_endstate_validation : bool, default True
                whether to conduct an endstate validation of the hybrid topology factory
-            ligand_file : str, default None
-                path to ligand of interest (i.e. small molecule or protein); .sdf or .pdb
+            ligand_input : str or list of oemol objects, default None


@glass-w list of oemol objects should be changed to oemol here

Ah yes good catch!

zhang-ivy · 2020-10-27T13:48:21Z

@glass-w have you re-run your example script above after making the changes today? just want to make sure that changing elif isinstance(ligand_input, oechem.OEMol): doesn't break anything

glass-w · 2020-10-27T16:21:43Z

@zhang-ivy yes I've re-run in a batch job using this functionality:

from perses.app.relative_point_mutation_setup import PointMutationExecutor
from perses.utils.openeye import generate_unique_atom_names
import oechem
from simtk import unit

ifs = oechem.oemolistream()
ifs.open('./ligand.sdf')
mol_list = [oechem.OEMol(mol) for mol in ifs.GetOEMols()]
mol = mol_list[0]

if len([atom.GetName() for atom in mol.GetAtoms()]) > len(set([atom.GetName() for atom in mol.GetAtoms()])):
        mol = generate_unique_atom_names(mol)
oechem.OEAssignAromaticFlags(mol, oechem.OEAroModelOpenEye)
oechem.OEAssignHybridization(mol)
oechem.OEAddExplicitHydrogens(mol)
oechem.OEPerceiveChiral(mol)

protein_file = './protein.pdb'

x = PointMutationExecutor(
protein_file,
'1',
'510', # an example residue number
'VAL', # an example mutation
ligand_input=mol,
ionic_strength=0.15 * unit.molar,
small_molecule_forcefields='openff-1.2.0'
)

it runs and produces all of the relvant output files. Again the above script is for my use case but for general OEMol preparation the user should use the lines you mentioned in your earlier comment.

dominicrufa · 2020-10-27T20:32:53Z

is there anything lossy between oemols and OFF molecule objects?

I believe some non-chemical information, such as SD tags, would be lost. But atom ordering and geometry should be preserved.

Is there any particular information you're worried about being lost, @dominicrufa? If so, we can check (or you can just ask on the OFF slack).

i think this is sufficient

dominicrufa

looks good to me

use list of oemol for PointMutationExecutor input

fca050a

zhang-ivy reviewed Oct 26, 2020

View reviewed changes

glass-w added 3 commits October 27, 2020 11:33

pass oemol object and remove list

30a05e6

update function description

06252df

remove redundant generate_unique_atom_names()

5f99606

glass-w requested a review from zhang-ivy October 27, 2020 12:22

zhang-ivy reviewed Oct 27, 2020

View reviewed changes

update function description

51868e1

Merge branch 'master' into fix-pm-lig

320eb7f

dominicrufa approved these changes Oct 27, 2020

View reviewed changes

zhang-ivy merged commit 7c5b932 into master Oct 27, 2020

zhang-ivy deleted the fix-pm-lig branch October 27, 2020 21:03

glass-w mentioned this pull request Oct 28, 2020

Refactor codebase to use Open Force Field Topology/Molecule objects #747

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Use a list of oemols for PointMutationExecutor input instead of files #746

Use a list of oemols for PointMutationExecutor input instead of files #746

glass-w commented Oct 26, 2020 •

edited

Loading

codecov bot commented Oct 26, 2020 •

edited

Loading

zhang-ivy left a comment •

edited

Loading

zhang-ivy Oct 26, 2020

dominicrufa Oct 26, 2020

zhang-ivy Oct 26, 2020

zhang-ivy Oct 26, 2020

zhang-ivy Oct 26, 2020

zhang-ivy Oct 26, 2020

zhang-ivy Oct 26, 2020

zhang-ivy Oct 26, 2020

dominicrufa Oct 26, 2020

zhang-ivy Oct 26, 2020

zhang-ivy commented Oct 26, 2020

jchodera commented Oct 26, 2020

dominicrufa commented Oct 26, 2020

jchodera commented Oct 27, 2020

glass-w commented Oct 27, 2020 •

edited

Loading

glass-w commented Oct 27, 2020

zhang-ivy Oct 27, 2020

glass-w Oct 27, 2020

zhang-ivy commented Oct 27, 2020

glass-w commented Oct 27, 2020 •

edited

Loading

dominicrufa commented Oct 27, 2020

dominicrufa left a comment

Use a list of oemols for PointMutationExecutor input instead of files #746

Use a list of oemols for PointMutationExecutor input instead of files #746

Conversation

glass-w commented Oct 26, 2020 • edited Loading

codecov bot commented Oct 26, 2020 • edited Loading

Codecov Report

zhang-ivy left a comment • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

zhang-ivy commented Oct 26, 2020

jchodera commented Oct 26, 2020

dominicrufa commented Oct 26, 2020

jchodera commented Oct 27, 2020

glass-w commented Oct 27, 2020 • edited Loading

glass-w commented Oct 27, 2020

Choose a reason for hiding this comment

Choose a reason for hiding this comment

zhang-ivy commented Oct 27, 2020

glass-w commented Oct 27, 2020 • edited Loading

dominicrufa commented Oct 27, 2020

dominicrufa left a comment

Choose a reason for hiding this comment

glass-w commented Oct 26, 2020 •

edited

Loading

codecov bot commented Oct 26, 2020 •

edited

Loading

zhang-ivy left a comment •

edited

Loading

glass-w commented Oct 27, 2020 •

edited

Loading

glass-w commented Oct 27, 2020 •

edited

Loading