Re: [Rdkit-discuss] Getting to grips with Open3DAlign

[Greg Landrum]

Hi Tim,

The below doesn't do what you think it does:

On Wed, Jun 22, 2016 at 10:29 AM, Tim Dudgeon  wrote:

This line:

AllChem.AddHs(mol1)
>

returns a new molecule with Hs added, it does not actually modify mol1.

The consequence of that is that this:

cids = AllChem.EmbedMultipleConfs(mol1, numConfs=500, maxAttempts=1000,
> pruneRmsThresh=0.1, useExpTorsionAnglePrefs=True, useBasicKnowledge=True,
> useRandomCoords=True, enforceChirality=True, numThreads=0)
>
>
Is still working on a molecule without Hs and will not give the best
conformations.

-greg
--
Attend Shape: An AT Tech Expo July 15-16. Meet us at AT Park in San
Francisco, CA to explore cutting-edge tech and listen to tech luminaries
present their vision of the future. This family event has something for
everyone, including kids. Get more information and register today.
http://sdm.link/attshape___
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss

Re: [Rdkit-discuss] Getting to grips with Open3DAlign

[Tim Dudgeon]

Paolo, Greg Sereina,

Thanks for your thoughts. In combination they do indeed resolve the 
problems.


To clarify one other thing: better O3A alignments give higher scores.

For the record here is my updated script:

from rdkit import Chem, RDConfig
from rdkit.Chem import AllChem, rdMolAlign

ref = 
Chem.MolFromSmiles('NC(=[NH2+])c1ccc(C[C@@H](NC(=O)CNS(=O)(=O)c2ccc3c3c2)C(=O)N2C2)cc1')
mol1 = 
Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1DWD_ligand.pdb')

mol1 = AllChem.AssignBondOrdersFromTemplate(ref, mol1)
mol2 = 
Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1PPC_ligand.pdb')

mol2 = AllChem.AssignBondOrdersFromTemplate(ref, mol2)

pyO3A = rdMolAlign.GetO3A(mol1, mol2)
pyO3A.Align()
score = pyO3A.Score()
print "Orig",score
Chem.MolToMolFile(mol1, "orig.mol")

AllChem.AddHs(mol1)
cids = AllChem.EmbedMultipleConfs(mol1, numConfs=500, maxAttempts=1000, 
pruneRmsThresh=0.1, useExpTorsionAnglePrefs=True, 
useBasicKnowledge=True, useRandomCoords=True, enforceChirality=True, 
numThreads=0)


pyO3As = rdMolAlign.GetO3AForProbeConfs(mol1, mol2, numThreads=0)
i = 0
highest = 0.0
for pyO3A in pyO3As:
pyO3A.Align()
score = pyO3A.Score()
if score > highest:
highest = score
highestConfId = i
i +=1

print "Highest:", highest, highestConfId
Chem.MolToMolFile(mol1, "highest.mol", confId=highestConfId)

I do have some other observations about the conformer generation, but 
I'll start a new topic for that.


Many thanks
Tim

On 22/06/2016 05:58, Sereina wrote:
I just tried Tim’s example (or a version of it that Greg sent me, 
respectively).
What is missing are the hydrogens for the torsion terms of ETKDG to 
work properly. Before generating the conformations AllChem.AddHs() 
should be called.


Best,
Sereina


On 22 Jun 2016, at 06:48, Sereina > wrote:


Based on the code snippets, Paolo has not used the basic-knowledge 
terms whereas Tim did.


When setting useExpTorsionAnglePrefs=True and useBasicKnowledge=True, 
a minimization is in principle not necessary anymore (unless there 
are aliphatic rings, because we currently don’t have torsion rules 
for them).


Best,
Sereina


On 22 Jun 2016, at 05:02, Greg Landrum > wrote:


I don't have anything to add to the pieces about the alignment 
(Paolo is the expert there!), but one comment on the conformation 
generation: If you used the background knowledge terms in the 
embedding, I don't think you should be getting the really distorted 
aromatic rings. Since in this case that does happen, for at least 
some conformations, I suspect there may be something wrong in the code.


I'll take a look at that (and ask Sereina too).

Best,
-greg


On Tue, Jun 21, 2016 at 10:30 PM, Paolo Tosco > wrote:


Dear Tim,

the Align() method returns an RMSD value, which however is
computed only on a limited number of atom pairs, namely those
that the algorithm was able to match between the two molecules,
so a low value is not particularly informative of the overall
goodness of the alignment, as it only indicates that the matched
atoms were matched nicely, but there might only be a few of
those in the core, while side chains are scattered all over.
The Score() method instead returns the O3AScore for the
alignment, which is a better way to assess the quality of the
superimposition.

The other problem in your script is that the i index is
incremented before recording it in the lowest/highest variables,
so the confIds are shifted by 1, as the conformation index in
the RDKit is 0-based.

I also noticed that without minimizing the conformations the
aromatic rings look quite distorted, so I added a MMFF
minimization, and I increased the number of generated
conformations and the pruneRmsThreshold. Setting to False the
experimental torsion angle preferences and basic knowledge about
rings seems to yield a larger variety of geometries which helps
reproducing this quite peculiar x-ray geometry which is probably
not so commonly found. Please find the modified script below.

Hope this helps, kind regards
Paolo


#!/usr/bin/env python


from rdkit import Chem, RDConfig
from rdkit.Chem import AllChem, rdMolAlign

ref =

Chem.MolFromSmiles('NC(=[NH2+])c1ccc(C[C@@H](NC(=O)CNS(=O)(=O)c2ccc3c3c2)C(=O)N2C2)cc1')
mol1 =

Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1DWD_ligand.pdb')
mol1 = AllChem.AssignBondOrdersFromTemplate(ref, mol1)
mol2 =

Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1PPC_ligand.pdb')
mol2 = AllChem.AssignBondOrdersFromTemplate(ref, mol2)

pyO3A = rdMolAlign.GetO3A(mol1, mol2)
rmsd = pyO3A.Align()
score = pyO3A.Score()
print "Orig",score
Chem.MolToMolFile(mol1, 

Re: [Rdkit-discuss] Getting to grips with Open3DAlign

[David Cosgrove]

Hi Sereina,

I beg to differ on the advisability of minimisation, even after using the
parameters you suggest to generate the conformation. I've recently been
using the CCDC's excellent Python API to analyse the results of the
generated conformations. This lets you very quickly assess whether any of
the bond distances or angles in a molecule is "unusual" with reference to
the data in the CSD. I've not done a systematic examination but my
anecdotal result is that you get far fewer such unusual geometries if you
give it a quick minimisation afterwards. Using the default optimiser
settings is quick, painless and, in my limited number of tests, worth the
effort.

I'd be interested to know if anyone else has different experience.

Cheers,
Dave


On Wednesday, 22 June 2016, Sereina  wrote:

> Based on the code snippets, Paolo has not used the basic-knowledge terms
> whereas Tim did.
>
> When setting useExpTorsionAnglePrefs=True and useBasicKnowledge=True, a
> minimization is in principle not necessary anymore (unless there are
> aliphatic rings, because we currently don’t have torsion rules for them).
>
> Best,
> Sereina
>
>
> On 22 Jun 2016, at 05:02, Greg Landrum  > wrote:
>
> I don't have anything to add to the pieces about the alignment (Paolo is
> the expert there!), but one comment on the conformation generation: If you
> used the background knowledge terms in the embedding, I don't think you
> should be getting the really distorted aromatic rings. Since in this case
> that does happen, for at least some conformations, I suspect there may be
> something wrong in the code.
>
> I'll take a look at that (and ask Sereina too).
>
> Best,
> -greg
>
>
> On Tue, Jun 21, 2016 at 10:30 PM, Paolo Tosco  > wrote:
>
>> Dear Tim,
>>
>> the Align() method returns an RMSD value, which however is computed only
>> on a limited number of atom pairs, namely those that the algorithm was able
>> to match between the two molecules, so a low value is not particularly
>> informative of the overall goodness of the alignment, as it only indicates
>> that the matched atoms were matched nicely, but there might only be a few
>> of those in the core, while side chains are scattered all over.
>> The Score() method instead returns the O3AScore for the alignment, which
>> is a better way to assess the quality of the superimposition.
>>
>> The other problem in your script is that the i index is incremented
>> before recording it in the lowest/highest variables, so the confIds are
>> shifted by 1, as the conformation index in the RDKit is 0-based.
>>
>> I also noticed that without minimizing the conformations the aromatic
>> rings look quite distorted, so I added a MMFF minimization, and I increased
>> the number of generated conformations and the pruneRmsThreshold. Setting to
>> False the experimental torsion angle preferences and basic knowledge about
>> rings seems to yield a larger variety of geometries which helps reproducing
>> this quite peculiar x-ray geometry which is probably not so commonly found.
>> Please find the modified script below.
>>
>> Hope this helps, kind regards
>> Paolo
>>
>>
>> #!/usr/bin/env python
>>
>>
>> from rdkit import Chem, RDConfig
>> from rdkit.Chem import AllChem, rdMolAlign
>>
>> ref = Chem.MolFromSmiles('NC(=[NH2+])c1ccc(C[C@
>> @H](NC(=O)CNS(=O)(=O)c2ccc3c3c2)C(=O)N2C2)cc1')
>> mol1 =
>> Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1DWD_ligand.pdb')
>> mol1 = AllChem.AssignBondOrdersFromTemplate(ref, mol1)
>> mol2 =
>> Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1PPC_ligand.pdb')
>> mol2 = AllChem.AssignBondOrdersFromTemplate(ref, mol2)
>>
>> pyO3A = rdMolAlign.GetO3A(mol1, mol2)
>> rmsd = pyO3A.Align()
>> score = pyO3A.Score()
>> print "Orig",score
>> Chem.MolToMolFile(mol1, "orig.mol")
>>
>> cids = AllChem.EmbedMultipleConfs(mol1, numConfs=250, maxAttempts=100,
>> pruneRmsThresh=0.5, useExpTorsionAnglePrefs=False,
>> useBasicKnowledge=False)
>> AllChem.MMFFOptimizeMoleculeConfs(mol1, mmffVariant='MMFF94s')
>> pyO3As = rdMolAlign.GetO3AForProbeConfs(mol1, mol2, numThreads=0)
>> i = 0
>> lowest = 9.9
>> highest = 0.0
>> for pyO3A in pyO3As:
>> rmsd = pyO3A.Align()
>> score = pyO3A.Score()
>> if score < lowest:
>> lowest = score
>> lowestConfId = i
>> if score > highest:
>> highest = score
>> highestConfId = i
>> i +=1
>>
>> print "Lowest:", lowest, lowestConfId
>> print "Highest:", highest, highestConfId
>>
>> Chem.MolToMolFile(mol1, "lowest.mol", confId=lowestConfId)
>> Chem.MolToMolFile(mol1, "highest.mol", confId=highestConfId)
>>
>>
>> On 06/21/16 15:41, Tim Dudgeon wrote:
>>
>> Hi All,
>>
>> I'm trying to get to grips with using Open3D Align in RDKit, but hitting
>> problems.
>>
>> My approach is to 

Re: [Rdkit-discuss] Getting to grips with Open3DAlign

[Sereina]

I just tried Tim’s example (or a version of it that Greg sent me, 
respectively). 
What is missing are the hydrogens for the torsion terms of ETKDG to work 
properly. Before generating the conformations AllChem.AddHs() should be called.

Best,
Sereina


On 22 Jun 2016, at 06:48, Sereina  wrote:

> Based on the code snippets, Paolo has not used the basic-knowledge terms 
> whereas Tim did. 
> 
> When setting useExpTorsionAnglePrefs=True and useBasicKnowledge=True, a 
> minimization is in principle not necessary anymore (unless there are 
> aliphatic rings, because we currently don’t have torsion rules for them).
> 
> Best,
> Sereina
> 
> 
> On 22 Jun 2016, at 05:02, Greg Landrum  wrote:
> 
>> I don't have anything to add to the pieces about the alignment (Paolo is the 
>> expert there!), but one comment on the conformation generation: If you used 
>> the background knowledge terms in the embedding, I don't think you should be 
>> getting the really distorted aromatic rings. Since in this case that does 
>> happen, for at least some conformations, I suspect there may be something 
>> wrong in the code.
>> 
>> I'll take a look at that (and ask Sereina too).
>> 
>> Best,
>> -greg
>> 
>> 
>> On Tue, Jun 21, 2016 at 10:30 PM, Paolo Tosco  wrote:
>> Dear Tim,
>> 
>> the Align() method returns an RMSD value, which however is computed only on 
>> a limited number of atom pairs, namely those that the algorithm was able to 
>> match between the two molecules, so a low value is not particularly 
>> informative of the overall goodness of the alignment, as it only indicates 
>> that the matched atoms were matched nicely, but there might only be a few of 
>> those in the core, while side chains are scattered all over.
>> The Score() method instead returns the O3AScore for the alignment, which is 
>> a better way to assess the quality of the superimposition.
>> 
>> The other problem in your script is that the i index is incremented before 
>> recording it in the lowest/highest variables, so the confIds are shifted by 
>> 1, as the conformation index in the RDKit is 0-based.
>> 
>> I also noticed that without minimizing the conformations the aromatic rings 
>> look quite distorted, so I added a MMFF minimization, and I increased the 
>> number of generated conformations and the pruneRmsThreshold. Setting to 
>> False the experimental torsion angle preferences and basic knowledge about 
>> rings seems to yield a larger variety of geometries which helps reproducing 
>> this quite peculiar x-ray geometry which is probably not so commonly found. 
>> Please find the modified script below.
>> 
>> Hope this helps, kind regards
>> Paolo
>> 
>> 
>> #!/usr/bin/env python
>> 
>> 
>> from rdkit import Chem, RDConfig
>> from rdkit.Chem import AllChem, rdMolAlign
>> 
>> ref = 
>> Chem.MolFromSmiles('NC(=[NH2+])c1ccc(C[C@@H](NC(=O)CNS(=O)(=O)c2ccc3c3c2)C(=O)N2C2)cc1')
>> mol1 = 
>> Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1DWD_ligand.pdb')
>> mol1 = AllChem.AssignBondOrdersFromTemplate(ref, mol1)
>> mol2 = 
>> Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1PPC_ligand.pdb')
>> mol2 = AllChem.AssignBondOrdersFromTemplate(ref, mol2)
>> 
>> pyO3A = rdMolAlign.GetO3A(mol1, mol2)
>> rmsd = pyO3A.Align()
>> score = pyO3A.Score()
>> print "Orig",score
>> Chem.MolToMolFile(mol1, "orig.mol")
>> 
>> cids = AllChem.EmbedMultipleConfs(mol1, numConfs=250, maxAttempts=100,
>> pruneRmsThresh=0.5, useExpTorsionAnglePrefs=False,
>> useBasicKnowledge=False)
>> AllChem.MMFFOptimizeMoleculeConfs(mol1, mmffVariant='MMFF94s')
>> pyO3As = rdMolAlign.GetO3AForProbeConfs(mol1, mol2, numThreads=0)
>> i = 0
>> lowest = 9.9
>> highest = 0.0
>> for pyO3A in pyO3As:
>> rmsd = pyO3A.Align()
>> score = pyO3A.Score()
>> if score < lowest:
>> lowest = score
>> lowestConfId = i
>> if score > highest:
>> highest = score
>> highestConfId = i
>> i +=1
>> 
>> print "Lowest:", lowest, lowestConfId
>> print "Highest:", highest, highestConfId
>> 
>> Chem.MolToMolFile(mol1, "lowest.mol", confId=lowestConfId)
>> Chem.MolToMolFile(mol1, "highest.mol", confId=highestConfId)
>> 
>> 
>> On 06/21/16 15:41, Tim Dudgeon wrote:
>>> Hi All,
>>> 
>>> I'm trying to get to grips with using Open3D Align in RDKit, but hitting 
>>> problems.
>>> 
>>> My approach is to generate random conformers of the probe molecule and 
>>> align it to the reference molecule.  My example is cobbled together from 
>>> the examples in the cookbook.
>>> 
>>> 
>>> 
>>> from rdkit import Chem, RDConfig
>>> from rdkit.Chem import AllChem, rdMolAlign
>>> 
>>> ref = 
>>> Chem.MolFromSmiles('NC(=[NH2+])c1ccc(C[C@@H](NC(=O)CNS(=O)(=O)c2ccc3c3c2)C(=O)N2C2)cc1')
>>> mol1 = 
>>> Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1DWD_ligand.pdb')
>>> mol1 = AllChem.AssignBondOrdersFromTemplate(ref, mol1)
>>> mol2 = 
>>> 

Re: [Rdkit-discuss] Getting to grips with Open3DAlign

[Sereina]

Based on the code snippets, Paolo has not used the basic-knowledge terms 
whereas Tim did. 

When setting useExpTorsionAnglePrefs=True and useBasicKnowledge=True, a 
minimization is in principle not necessary anymore (unless there are aliphatic 
rings, because we currently don’t have torsion rules for them).

Best,
Sereina


On 22 Jun 2016, at 05:02, Greg Landrum  wrote:

> I don't have anything to add to the pieces about the alignment (Paolo is the 
> expert there!), but one comment on the conformation generation: If you used 
> the background knowledge terms in the embedding, I don't think you should be 
> getting the really distorted aromatic rings. Since in this case that does 
> happen, for at least some conformations, I suspect there may be something 
> wrong in the code.
> 
> I'll take a look at that (and ask Sereina too).
> 
> Best,
> -greg
> 
> 
> On Tue, Jun 21, 2016 at 10:30 PM, Paolo Tosco  wrote:
> Dear Tim,
> 
> the Align() method returns an RMSD value, which however is computed only on a 
> limited number of atom pairs, namely those that the algorithm was able to 
> match between the two molecules, so a low value is not particularly 
> informative of the overall goodness of the alignment, as it only indicates 
> that the matched atoms were matched nicely, but there might only be a few of 
> those in the core, while side chains are scattered all over.
> The Score() method instead returns the O3AScore for the alignment, which is a 
> better way to assess the quality of the superimposition.
> 
> The other problem in your script is that the i index is incremented before 
> recording it in the lowest/highest variables, so the confIds are shifted by 
> 1, as the conformation index in the RDKit is 0-based.
> 
> I also noticed that without minimizing the conformations the aromatic rings 
> look quite distorted, so I added a MMFF minimization, and I increased the 
> number of generated conformations and the pruneRmsThreshold. Setting to False 
> the experimental torsion angle preferences and basic knowledge about rings 
> seems to yield a larger variety of geometries which helps reproducing this 
> quite peculiar x-ray geometry which is probably not so commonly found. Please 
> find the modified script below.
> 
> Hope this helps, kind regards
> Paolo
> 
> 
> #!/usr/bin/env python
> 
> 
> from rdkit import Chem, RDConfig
> from rdkit.Chem import AllChem, rdMolAlign
> 
> ref = 
> Chem.MolFromSmiles('NC(=[NH2+])c1ccc(C[C@@H](NC(=O)CNS(=O)(=O)c2ccc3c3c2)C(=O)N2C2)cc1')
> mol1 = 
> Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1DWD_ligand.pdb')
> mol1 = AllChem.AssignBondOrdersFromTemplate(ref, mol1)
> mol2 = 
> Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1PPC_ligand.pdb')
> mol2 = AllChem.AssignBondOrdersFromTemplate(ref, mol2)
> 
> pyO3A = rdMolAlign.GetO3A(mol1, mol2)
> rmsd = pyO3A.Align()
> score = pyO3A.Score()
> print "Orig",score
> Chem.MolToMolFile(mol1, "orig.mol")
> 
> cids = AllChem.EmbedMultipleConfs(mol1, numConfs=250, maxAttempts=100,
> pruneRmsThresh=0.5, useExpTorsionAnglePrefs=False,
> useBasicKnowledge=False)
> AllChem.MMFFOptimizeMoleculeConfs(mol1, mmffVariant='MMFF94s')
> pyO3As = rdMolAlign.GetO3AForProbeConfs(mol1, mol2, numThreads=0)
> i = 0
> lowest = 9.9
> highest = 0.0
> for pyO3A in pyO3As:
> rmsd = pyO3A.Align()
> score = pyO3A.Score()
> if score < lowest:
> lowest = score
> lowestConfId = i
> if score > highest:
> highest = score
> highestConfId = i
> i +=1
> 
> print "Lowest:", lowest, lowestConfId
> print "Highest:", highest, highestConfId
> 
> Chem.MolToMolFile(mol1, "lowest.mol", confId=lowestConfId)
> Chem.MolToMolFile(mol1, "highest.mol", confId=highestConfId)
> 
> 
> On 06/21/16 15:41, Tim Dudgeon wrote:
>> Hi All,
>> 
>> I'm trying to get to grips with using Open3D Align in RDKit, but hitting 
>> problems.
>> 
>> My approach is to generate random conformers of the probe molecule and align 
>> it to the reference molecule.  My example is cobbled together from the 
>> examples in the cookbook.
>> 
>> 
>> 
>> from rdkit import Chem, RDConfig
>> from rdkit.Chem import AllChem, rdMolAlign
>> 
>> ref = 
>> Chem.MolFromSmiles('NC(=[NH2+])c1ccc(C[C@@H](NC(=O)CNS(=O)(=O)c2ccc3c3c2)C(=O)N2C2)cc1')
>> mol1 = 
>> Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1DWD_ligand.pdb')
>> mol1 = AllChem.AssignBondOrdersFromTemplate(ref, mol1)
>> mol2 = 
>> Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1PPC_ligand.pdb')
>> mol2 = AllChem.AssignBondOrdersFromTemplate(ref, mol2)
>> 
>> pyO3A = rdMolAlign.GetO3A(mol1, mol2)
>> score = pyO3A.Align()
>> print "Orig",score
>> Chem.MolToMolFile(mol1, "orig.mol")
>> 
>> cids = AllChem.EmbedMultipleConfs(mol1, numConfs=100, maxAttempts=100, 
>> pruneRmsThresh=0.1, useExpTorsionAnglePrefs=True, useBasicKnowledge=True)
>> 
>> pyO3As = 

Re: [Rdkit-discuss] Getting to grips with Open3DAlign

[Greg Landrum]

I don't have anything to add to the pieces about the alignment (Paolo is
the expert there!), but one comment on the conformation generation: If you
used the background knowledge terms in the embedding, I don't think you
should be getting the really distorted aromatic rings. Since in this case
that does happen, for at least some conformations, I suspect there may be
something wrong in the code.

I'll take a look at that (and ask Sereina too).

Best,
-greg


On Tue, Jun 21, 2016 at 10:30 PM, Paolo Tosco  wrote:

> Dear Tim,
>
> the Align() method returns an RMSD value, which however is computed only
> on a limited number of atom pairs, namely those that the algorithm was able
> to match between the two molecules, so a low value is not particularly
> informative of the overall goodness of the alignment, as it only indicates
> that the matched atoms were matched nicely, but there might only be a few
> of those in the core, while side chains are scattered all over.
> The Score() method instead returns the O3AScore for the alignment, which
> is a better way to assess the quality of the superimposition.
>
> The other problem in your script is that the i index is incremented before
> recording it in the lowest/highest variables, so the confIds are shifted by
> 1, as the conformation index in the RDKit is 0-based.
>
> I also noticed that without minimizing the conformations the aromatic
> rings look quite distorted, so I added a MMFF minimization, and I increased
> the number of generated conformations and the pruneRmsThreshold. Setting to
> False the experimental torsion angle preferences and basic knowledge about
> rings seems to yield a larger variety of geometries which helps reproducing
> this quite peculiar x-ray geometry which is probably not so commonly found.
> Please find the modified script below.
>
> Hope this helps, kind regards
> Paolo
>
>
> #!/usr/bin/env python
>
>
> from rdkit import Chem, RDConfig
> from rdkit.Chem import AllChem, rdMolAlign
>
> ref = Chem.MolFromSmiles('NC(=[NH2+])c1ccc(C[C@
> @H](NC(=O)CNS(=O)(=O)c2ccc3c3c2)C(=O)N2C2)cc1')
> mol1 =
> Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1DWD_ligand.pdb')
> mol1 = AllChem.AssignBondOrdersFromTemplate(ref, mol1)
> mol2 =
> Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1PPC_ligand.pdb')
> mol2 = AllChem.AssignBondOrdersFromTemplate(ref, mol2)
>
> pyO3A = rdMolAlign.GetO3A(mol1, mol2)
> rmsd = pyO3A.Align()
> score = pyO3A.Score()
> print "Orig",score
> Chem.MolToMolFile(mol1, "orig.mol")
>
> cids = AllChem.EmbedMultipleConfs(mol1, numConfs=250, maxAttempts=100,
> pruneRmsThresh=0.5, useExpTorsionAnglePrefs=False,
> useBasicKnowledge=False)
> AllChem.MMFFOptimizeMoleculeConfs(mol1, mmffVariant='MMFF94s')
> pyO3As = rdMolAlign.GetO3AForProbeConfs(mol1, mol2, numThreads=0)
> i = 0
> lowest = 9.9
> highest = 0.0
> for pyO3A in pyO3As:
> rmsd = pyO3A.Align()
> score = pyO3A.Score()
> if score < lowest:
> lowest = score
> lowestConfId = i
> if score > highest:
> highest = score
> highestConfId = i
> i +=1
>
> print "Lowest:", lowest, lowestConfId
> print "Highest:", highest, highestConfId
>
> Chem.MolToMolFile(mol1, "lowest.mol", confId=lowestConfId)
> Chem.MolToMolFile(mol1, "highest.mol", confId=highestConfId)
>
>
> On 06/21/16 15:41, Tim Dudgeon wrote:
>
> Hi All,
>
> I'm trying to get to grips with using Open3D Align in RDKit, but hitting
> problems.
>
> My approach is to generate random conformers of the probe molecule and
> align it to the reference molecule.  My example is cobbled together from
> the examples in the cookbook.
>
> from rdkit import Chem, RDConfig
> from rdkit.Chem import AllChem, rdMolAlign
>
> ref = Chem.MolFromSmiles('NC(=[NH2+])c1ccc(C[C@
> @H](NC(=O)CNS(=O)(=O)c2ccc3c3c2)C(=O)N2C2)cc1')
> mol1 =
> Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1DWD_ligand.pdb')
> mol1 = AllChem.AssignBondOrdersFromTemplate(ref, mol1)
> mol2 =
> Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1PPC_ligand.pdb')
> mol2 = AllChem.AssignBondOrdersFromTemplate(ref, mol2)
>
> pyO3A = rdMolAlign.GetO3A(mol1, mol2)
> score = pyO3A.Align()
> print "Orig",score
> Chem.MolToMolFile(mol1, "orig.mol")
>
> cids = AllChem.EmbedMultipleConfs(mol1, numConfs=100, maxAttempts=100,
> pruneRmsThresh=0.1, useExpTorsionAnglePrefs=True, useBasicKnowledge=True)
>
> pyO3As = rdMolAlign.GetO3AForProbeConfs(mol1, mol2, numThreads=0)
> i = 0
> lowest = 9.9
> highest = 0.0
> for pyO3A in pyO3As:
> i +=1
> score = pyO3A.Align()
> if score < lowest:
> lowest = score
> lowestConfId = i
> if score > highest:
> highest = score
> highestConfId = i
>
> print "Lowest:", lowest, lowestConfId
> print "Highest:", highest, highestConfId
>
> Chem.MolToMolFile(mol1, "lowest.mol", confId=lowestConfId)
> Chem.MolToMolFile(mol1, "highest.mol", 

Re: [Rdkit-discuss] Getting to grips with Open3DAlign

[Paolo Tosco]

Dear Tim,

the Align() method returns an RMSD value, which however is computed only 
on a limited number of atom pairs, namely those that the algorithm was 
able to match between the two molecules, so a low value is not 
particularly informative of the overall goodness of the alignment, as it 
only indicates that the matched atoms were matched nicely, but there 
might only be a few of those in the core, while side chains are 
scattered all over.
The Score() method instead returns the O3AScore for the alignment, which 
is a better way to assess the quality of the superimposition.


The other problem in your script is that the i index is incremented 
before recording it in the lowest/highest variables, so the confIds are 
shifted by 1, as the conformation index in the RDKit is 0-based.


I also noticed that without minimizing the conformations the aromatic 
rings look quite distorted, so I added a MMFF minimization, and I 
increased the number of generated conformations and the 
pruneRmsThreshold. Setting to False the experimental torsion angle 
preferences and basic knowledge about rings seems to yield a larger 
variety of geometries which helps reproducing this quite peculiar x-ray 
geometry which is probably not so commonly found. Please find the 
modified script below.


Hope this helps, kind regards
Paolo


#!/usr/bin/env python


from rdkit import Chem, RDConfig
from rdkit.Chem import AllChem, rdMolAlign

ref = 
Chem.MolFromSmiles('NC(=[NH2+])c1ccc(C[C@@H](NC(=O)CNS(=O)(=O)c2ccc3c3c2)C(=O)N2C2)cc1')
mol1 = 
Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1DWD_ligand.pdb')

mol1 = AllChem.AssignBondOrdersFromTemplate(ref, mol1)
mol2 = 
Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1PPC_ligand.pdb')

mol2 = AllChem.AssignBondOrdersFromTemplate(ref, mol2)

pyO3A = rdMolAlign.GetO3A(mol1, mol2)
rmsd = pyO3A.Align()
score = pyO3A.Score()
print "Orig",score
Chem.MolToMolFile(mol1, "orig.mol")

cids = AllChem.EmbedMultipleConfs(mol1, numConfs=250, maxAttempts=100,
pruneRmsThresh=0.5, useExpTorsionAnglePrefs=False,
useBasicKnowledge=False)
AllChem.MMFFOptimizeMoleculeConfs(mol1, mmffVariant='MMFF94s')
pyO3As = rdMolAlign.GetO3AForProbeConfs(mol1, mol2, numThreads=0)
i = 0
lowest = 9.9
highest = 0.0
for pyO3A in pyO3As:
rmsd = pyO3A.Align()
score = pyO3A.Score()
if score < lowest:
lowest = score
lowestConfId = i
if score > highest:
highest = score
highestConfId = i
i +=1

print "Lowest:", lowest, lowestConfId
print "Highest:", highest, highestConfId

Chem.MolToMolFile(mol1, "lowest.mol", confId=lowestConfId)
Chem.MolToMolFile(mol1, "highest.mol", confId=highestConfId)


On 06/21/16 15:41, Tim Dudgeon wrote:


Hi All,

I'm trying to get to grips with using Open3D Align in RDKit, but 
hitting problems.


My approach is to generate random conformers of the probe molecule and 
align it to the reference molecule.  My example is cobbled together 
from the examples in the cookbook.


from rdkit import Chem, RDConfig
from rdkit.Chem import AllChem, rdMolAlign

ref = 
Chem.MolFromSmiles('NC(=[NH2+])c1ccc(C[C@@H](NC(=O)CNS(=O)(=O)c2ccc3c3c2)C(=O)N2C2)cc1')
mol1 = 
Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1DWD_ligand.pdb')

mol1 = AllChem.AssignBondOrdersFromTemplate(ref, mol1)
mol2 = 
Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1PPC_ligand.pdb')

mol2 = AllChem.AssignBondOrdersFromTemplate(ref, mol2)

pyO3A = rdMolAlign.GetO3A(mol1, mol2)
score = pyO3A.Align()
print "Orig",score
Chem.MolToMolFile(mol1, "orig.mol")

cids = AllChem.EmbedMultipleConfs(mol1, numConfs=100, maxAttempts=100, 
pruneRmsThresh=0.1, useExpTorsionAnglePrefs=True, useBasicKnowledge=True)


pyO3As = rdMolAlign.GetO3AForProbeConfs(mol1, mol2, numThreads=0)
i = 0
lowest = 9.9
highest = 0.0
for pyO3A in pyO3As:
i +=1
score = pyO3A.Align()
if score < lowest:
lowest = score
lowestConfId = i
if score > highest:
highest = score
highestConfId = i

print "Lowest:", lowest, lowestConfId
print "Highest:", highest, highestConfId

Chem.MolToMolFile(mol1, "lowest.mol", confId=lowestConfId)
Chem.MolToMolFile(mol1, "highest.mol", confId=highestConfId)

What I'm finding is that the alignments with the lowest and highest 
O3A scores are much worse alignments (visually) than the original 
structure (the structure from 1DWD). Typical scores are:


Original 1DWD structure: 0.38
Lowest scoring conformer: 0.186
Highest scoring conformer: 0.78

Now I'm assuming that lower O3A align scores are better (though that's 
not specifically stated), so as my lowest score is lower than the 
original alignment I would have expected it to be a better alignment, 
but its clearly much worse. And if I'm wrong and higher scores are 
better then the same applies.


Clearly I've not understood something correctly!

Tim

[Rdkit-discuss] Getting to grips with Open3DAlign

[Tim Dudgeon]

Hi All,

I'm trying to get to grips with using Open3D Align in RDKit, but hitting 
problems.


My approach is to generate random conformers of the probe molecule and 
align it to the reference molecule.  My example is cobbled together from 
the examples in the cookbook.


from rdkit import Chem, RDConfig
from rdkit.Chem import AllChem, rdMolAlign

ref = 
Chem.MolFromSmiles('NC(=[NH2+])c1ccc(C[C@@H](NC(=O)CNS(=O)(=O)c2ccc3c3c2)C(=O)N2C2)cc1')
mol1 = 
Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1DWD_ligand.pdb')

mol1 = AllChem.AssignBondOrdersFromTemplate(ref, mol1)
mol2 = 
Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1PPC_ligand.pdb')

mol2 = AllChem.AssignBondOrdersFromTemplate(ref, mol2)

pyO3A = rdMolAlign.GetO3A(mol1, mol2)
score = pyO3A.Align()
print "Orig",score
Chem.MolToMolFile(mol1, "orig.mol")

cids = AllChem.EmbedMultipleConfs(mol1, numConfs=100, maxAttempts=100, 
pruneRmsThresh=0.1, useExpTorsionAnglePrefs=True, useBasicKnowledge=True)


pyO3As = rdMolAlign.GetO3AForProbeConfs(mol1, mol2, numThreads=0)
i = 0
lowest = 9.9
highest = 0.0
for pyO3A in pyO3As:
i +=1
score = pyO3A.Align()
if score < lowest:
lowest = score
lowestConfId = i
if score > highest:
highest = score
highestConfId = i

print "Lowest:", lowest, lowestConfId
print "Highest:", highest, highestConfId

Chem.MolToMolFile(mol1, "lowest.mol", confId=lowestConfId)
Chem.MolToMolFile(mol1, "highest.mol", confId=highestConfId)

What I'm finding is that the alignments with the lowest and highest O3A 
scores are much worse alignments (visually) than the original structure 
(the structure from 1DWD). Typical scores are:


Original 1DWD structure: 0.38
Lowest scoring conformer: 0.186
Highest scoring conformer: 0.78

Now I'm assuming that lower O3A align scores are better (though that's 
not specifically stated), so as my lowest score is lower than the 
original alignment I would have expected it to be a better alignment, 
but its clearly much worse. And if I'm wrong and higher scores are 
better then the same applies.


Clearly I've not understood something correctly!

Tim

--
Attend Shape: An AT Tech Expo July 15-16. Meet us at AT Park in San
Francisco, CA to explore cutting-edge tech and listen to tech luminaries
present their vision of the future. This family event has something for
everyone, including kids. Get more information and register today.
http://sdm.link/attshape___
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss