Dear Stéphane,
the reason why the rmsd is not computed for some of the structures is
that the chirality of atom 15 is inverted in carbLactosamineTypeI
compared to carbLactosamine, and since in the substructure match you
request useChirality=True, the stereochemistry doesn't match the PDB
reference anymore (see attached script).
If you either use useChirality=False, or you epimerize atom 15
uncommenting the SetChiralTag() command in the script, the rmsd will be
computed.
Cheers,
p.
On 11/17/2016 07:42 PM, Téletchéa Stéphane wrote:
Le 17/11/2016 à 20:12, Paolo Tosco a écrit :
Dear Stéphane,
I'll look into this tonight and let you know.
Best,
p.
Thanks a lot.
To add more input, depending on the input file provide, the rmsd is
computed or not ...
(change Chem.MolFromPDBFile(sys.argv[1]) for
Chem.MolFromMolFile(sys.argv[1]))
./calc_rmsd_simple.py sd_out/carbLactosamineTypeI5.sdf jeu_test.sd
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamine 0.000000
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamine 0.000000
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamine 0.000000
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamine 0.000000
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamine 0.000000
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamine 0.000000
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamine 0.000000
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamine 0.000000
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamine 0.000000
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamine 0.000000
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamineTypeI
0.669621
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamineTypeI
0.766064
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamineTypeI
0.912119
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamineTypeI
0.120689
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamineTypeI
0.512118
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamineTypeI
0.000000
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamineTypeI
1.899627
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamineTypeI
1.687619
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamineTypeI
0.000000
sd_out/carbLactosamineTypeI5.sdf jeu_test.sd carbLactosamineTypeI
0.000000
./calc_rmsd_simple.py sd_out/carbLactosamine5.sdf jeu_test.sd
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamine 0.680751
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamine 0.264023
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamine 0.127012
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamine 0.196644
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamine 0.176413
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamine 0.000000
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamine 0.578994
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamine 0.567569
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamine 0.592443
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamine 1.295632
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamineTypeI 0.000000
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamineTypeI 0.000000
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamineTypeI 0.000000
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamineTypeI 0.000000
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamineTypeI 0.000000
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamineTypeI 0.000000
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamineTypeI 0.000000
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamineTypeI 0.000000
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamineTypeI 0.000000
sd_out/carbLactosamine5.sdf jeu_test.sd carbLactosamineTypeI 0.000000
Very weird ... so this is not only an orientation issue ...
Stéphane
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import re
import sys
import os
import numpy as np
from collections import defaultdict
from rdkit import Chem
from rdkit.Chem import AllChem, SanitizeMol
from rdkit.Chem import SDMolSupplier, SDWriter
from rdkit.Chem import MolFromMol2File, RemoveHs
def getCoords(ligand):
coords=[]
for i in range(0, ligand.GetNumAtoms()):
pos = ligand.GetConformer().GetAtomPosition(i)
coords.append(pos)
return coords
def rmsd(coord1, coord2):
r=0
try:
r=np.sqrt( np.sum((np.array(coord1)-np.array(coord2))**2)/len(coord1))
except:
pass
return r
def main():
ref_gal = Chem.MolFromPDBFile(sys.argv[1])
Chem.AssignAtomChiralTagsFromStructure(ref_gal)
galcoords=getCoords(ref_gal)
try:
poses=open(sys.argv[2],'r')
except:
print "Unable to read %s" % sys.argv[2]
exit
suppl = SDMolSupplier(sys.argv[2])
i=0
chiralDict={
Chem.CHI_UNSPECIFIED: 0,
Chem.CHI_TETRAHEDRAL_CW: 1,
Chem.CHI_TETRAHEDRAL_CCW: -1
}
for l in suppl:
Chem.AssignAtomChiralTagsFromStructure(l,replaceExistingTags=True)
#l.GetAtomWithIdx(15).SetChiralTag(Chem.CHI_TETRAHEDRAL_CCW)
i=i+1
nom='mol'+str(i)
try:
nom=l.GetProp('_Name')
except:
pass
# first get all the atom coordinates from substructure match using chirality otherwise it won't work
m=l.GetSubstructMatches(ref_gal,useChirality=True)
#print m
l_coords=[]
if m:
for a in m[0]:
l_coords.append(l.GetConformer().GetAtomPosition(a))
# Second compute the rmsd to the reference ligand
r=rmsd(galcoords,l_coords)
try:
c = []
for a in l.GetAtoms():
tag = chiralDict[a.GetChiralTag()]
if (tag):
c.append((a.GetIdx(), tag))
print "%s %s %s %f %s %s" % (sys.argv[1], sys.argv[2], nom, r, str(m), str(c))
except Exception as e:
print "%s: unable to compute the rmsd (error: %s)" % (nom,e)
pass
if __name__ == "__main__":
main()
------------------------------------------------------------------------------
_______________________________________________
Rdkit-discuss mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
