Dear James,
On Fri, Jul 2, 2010 at 3:05 PM, James Davidson <j.david...@vernalis.com> wrote:
> Dear All,
>
> I am trying to work out the best way to accomplish some tasks involving
> RDKit, using PyMOL as an interface, and would appreciate some help. I would
> like to be able to start from a PDB file of a ligand-bound crystal structure
> loaded in PyMOL and be able to 'virtually' build some analogues - initially
> just simple substitutions - and visually inspect the results.
>
> (1) So my first question is - having started PyMOL with the -R option, is
> there an easy or recommended way of transferring molecules from PyMOL to
> RDKit? I can accomplish this by writing molfiles to a temporary file, but
> wondered if I am creating work, if eg RDKit could automatically create Mol
> objects from non-biopolymer atoms in PyMol(?). ie it would be nice if:
>
> from rdkit import Chem
> from rdkit.Chem import PyMol
> v = PyMol.MolViewer()
>
> # Invented function to create an RDKit mol object
> mol = v.GetAtomsAsMol(selection)
The only way I currently know of to get molecule data out of pymol is
to use the save command and a temporary file; one could create the
function above on the basis of the fact that this command works
without problems:
v.server.do('save /tmp/blah.mol,(sele)')
If there were sufficient demand for it, I'd guess that Jason would add
(or arrange to have added) an option to the save command to make the
results available as text instead of sending them to a file. It might
be worth asking him.
> (2) Once the ligand is converted to an RDKit mol object (by whatever means)
> I want to enumrate some libraries of virtual products - eg choosing an atom
> in PyMol as the attachment point, then running a set of reactions to get
> products with a set of substituents added. In principle I think this is
> quite straightforward; however, what I am struggling with is a mechanism to
> 'freeze' the 3D coordintes of the original ligand atoms, but still be able
> to use RDKit to generate sensible 3D positions for the newly added atoms so
> that the products can be passed back to PyMOL and minimised in situ (if
> required) using 'mengine.exe'. Am I looking at this the wrong way, and
> should I actually try aligning the virtual products back on the starting
> ligand conformation?
Well, you may be looking at it the wrong way, but I've looked at it
that way too, so we're at least wrong together. ;-)
There is an underdocumented (where, in this case under == un) function
that does part of what I believe you want in AllChem:
#---------------------------
# start by getting a geometry for our core, for the purposes of the demo:
[28]>>> core = Chem.MolFromSmiles('c1cc(C)cc2c1[nH]cc2')
[29]>>> core = Chem.AddHs(core)
[30]>>> AllChem.EmbedMolecule(core)
Out[30] 0
[31]>>> AllChem.UFFOptimizeMolecule(core)
Out[31] 0
# remove the Hs, otherwise we won't get a substruct match:
[32]>>> core = Chem.RemoveHs(core)
#---------------------------
# now generate the query molecule, with Hs :
[34]>>> mol = Chem.MolFromSmiles('CCc1cc(C)cc2c1[n](C(=O)OC)cc2')
[35]>>> mol = Chem.AddHs(mol)
# do the constrained embedding:
[37]>>> newMol=AllChem.ConstrainedEmbed(mol,core,True)
# the molecule has a property with the RMS for the final alignment:
[38]>>> mol.GetProp('EmbedRMS')
Out[38] '0.0278532402393'
# but we can also demonstrate that the geometries are similar by inspection:
[39]>>> print Chem.MolToMolBlock(core)
------> print(Chem.MolToMolBlock(core))
RDKit 3D
10 11 0 0 0 0 0 0 0 0999 V2000
-0.2038 1.6188 -0.4456 C 0 0 0 0 0 0 0 0 0 0 0 0
1.1538 1.2875 -0.3158 C 0 0 0 0 0 0 0 0 0 0 0 0
1.5492 -0.0245 0.0270 C 0 0 0 0 0 0 0 0 0 0 0 0
3.0012 -0.3502 0.2035 C 0 0 0 0 0 0 0 0 0 0 0 0
0.5703 -1.0078 0.2731 C 0 0 0 0 0 0 0 0 0 0 0 0
-0.7746 -0.6593 0.1391 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.1507 0.6210 -0.2130 C 0 0 0 0 0 0 0 0 0 0 0 0
-2.5018 0.7120 -0.2748 N 0 0 0 0 0 0 0 0 0 0 0 0
-2.9792 -0.5229 0.0418 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.9314 -1.3926 0.3038 C 0 0 0 0 0 0 0 0 0 0 0 0
1 2 2 0
2 3 1 0
3 4 1 0
3 5 2 0
5 6 1 0
6 7 2 0
7 8 1 0
8 9 1 0
9 10 2 0
7 1 1 0
10 6 1 0
M END
[40]>>> print Chem.MolToMolBlock(Chem.RemoveHs(newMol))
------> print(Chem.MolToMolBlock(Chem.RemoveHs(newMol)))
RDKit 3D
16 17 0 0 0 0 0 0 0 0999 V2000
0.0460 3.9887 0.5472 C 0 0 0 0 0 0 0 0 0 0 0 0
-0.4815 3.1257 -0.5989 C 0 0 0 0 0 0 0 0 0 0 0 0
-0.1879 1.6533 -0.4169 C 0 0 0 0 0 0 0 0 0 0 0 0
1.1668 1.2744 -0.3174 C 0 0 0 0 0 0 0 0 0 0 0 0
1.5491 -0.0364 0.0192 C 0 0 0 0 0 0 0 0 0 0 0 0
3.0001 -0.3651 0.2010 C 0 0 0 0 0 0 0 0 0 0 0 0
0.5630 -1.0025 0.2738 C 0 0 0 0 0 0 0 0 0 0 0 0
-0.7759 -0.6379 0.1381 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.1627 0.6457 -0.2263 C 0 0 0 0 0 0 0 0 0 0 0 0
-2.5280 0.6887 -0.3026 N 0 0 0 0 0 0 0 0 0 0 0 0
-3.4310 1.6754 -0.8441 C 0 0 0 0 0 0 0 0 0 0 0 0
-3.1599 2.2139 -2.0006 O 0 0 0 0 0 0 0 0 0 0 0 0
-4.7497 1.7964 -0.2732 O 0 0 0 0 0 0 0 0 0 0 0 0
-5.6798 2.7222 -0.9267 C 0 0 0 0 0 0 0 0 0 0 0 0
-2.9758 -0.5480 0.0511 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.9156 -1.3901 0.3193 C 0 0 0 0 0 0 0 0 0 0 0 0
1 2 1 0
2 3 1 0
3 4 2 0
4 5 1 0
5 6 1 0
5 7 2 0
7 8 1 0
8 9 2 0
9 10 1 0
11 12 2 3
11 13 1 0
13 14 1 0
10 11 1 0
10 15 1 0
15 16 2 0
9 3 1 0
16 8 1 0
M END
Of course these are distance-geometry conformations that have been
optimized with UFF, so they aren't the nicest things the world has
ever seen, but they do provide starting points for a real force field.
> (3) Apologies that this last point is maybe a bit off-topic, but I wondered
> if anyone has an opinion as to whether MMTK is the way to go for 'simple'
> minimisations of modified ligands bound to proteins? (I don't have any real
> experience with MMTK...
can't really help here, sorry
-greg
------------------------------------------------------------------------------
This SF.net email is sponsored by Sprint
What will you do first with EVO, the first 4G phone?
Visit sprint.com/first -- http://p.sf.net/sfu/sprint-com-first
_______________________________________________
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
