7
$\begingroup$

I recently started using both pysmiles and RDkit to parse SMILES strings into molecules. However, I sometimes got different results between the two libraries.

For example, on the molecule described by the string OCCn2c(=N)n(CCOc1ccc(Cl)cc1Cl)c3ccccc23, which is parsed using RDkit into the following molecule:

SMILES string converted to molecule by RDkit

This molecule has some atom such that RDkit's atom.GetNumImplicitHs() gives a value of $0$, while pysmiles hcount attribute for these atoms is $1$.

This also happens for other molecules, and also with the formal charge value I have seen a similar problem.

Why does this happen? How can I know which library is "right"? Or maybe I didn't understand correctly those features (I'm no chemist, so excuse me for that)? Thanks!

For completeness, the code I used to process the SMILES string with each library is shown below:

  • RDkit: molecule = RDkit.Chem.MolFromSmiles(smiles)
  • pysmiles: molecule = pysmiles.read_smiles(smiles)
$\endgroup$
6
  • $\begingroup$ How are you reading in the SMILES string with each program? They will classify hydrogens as implicit/explicit depending on the options passed to these functions. $\endgroup$
    – Tyberius
    Jul 29 at 0:19
  • $\begingroup$ Using RDkit, this is the code: molecule = RDkit.Chem.MolFromSmiles(smiles), and for pysmiles this is what I use: molecule = pysmiles.read_smiles(smiles). I don't have any additional options that I passed to them $\endgroup$
    – nir shahar
    Jul 29 at 7:19
  • $\begingroup$ is the pysmiles you are talking about the graph-to-smiles program on github? If so i’m wondering why you are using this like this. $\endgroup$
    – Cody Aldaz
    Jul 29 at 14:32
  • $\begingroup$ This is the library im using: pypi.org/project/pysmiles $\endgroup$
    – nir shahar
    Jul 29 at 14:34
  • $\begingroup$ Also, I'm converting from smiles to a graph, and not vice-versa $\endgroup$
    – nir shahar
    Jul 29 at 14:49
5
$\begingroup$

It is a bug in pysmiles.

pysmiles has problems with aromatic nitrogen if the SMILES is in aromatic form (lower case) but not if it is kekulized (upper case and double bonds)

If we add the hydrogens to the molecule, the nitrogens on index 3 and 6 have none, because they have already three bonds.

from rdkit import Chem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
import pysmiles

s1 = 'OCCn2c(=N)n(CCOc1ccc(Cl)cc1Cl)c3ccccc23'  # aromatic
s2 = 'OCCN2C(=N)N(CCOC1=CC=C(Cl)C=C1Cl)C3=CC=CC=C23'  # kekulized


def show_implicit_h(smiles):
    m = Chem.MolFromSmiles(smiles)
    for atom in m.GetAtoms():
        atom.SetProp('atomLabel', str(atom.GetIdx()))
    m = Chem.AddHs(m)
    return Draw.MolToImage(m, size=(300, 300))


show_implicit_h(s1)

aromatic

show_implicit_h(s2)

kekulized

Counting the hydrogen on both smiles shows that pysmiles finds hydrogen in the aromatic, but not in the kekulized SMILES. RDKit finds none in both.

Columns: Index, Symbol, hcount pysmiles, hcount RDKit

def count_implicit_h(s):
    m = Chem.MolFromSmiles(s)

    mol = pysmiles.read_smiles(s)
    hc = mol.nodes(data='hcount')

    for n in hc:
        print(n[0],
              m.GetAtomWithIdx(n[0]).GetSymbol(), n[1],
              m.GetAtomWithIdx(n[0]).GetNumImplicitHs())

aromatic form (count_implicit_h(s1)):

0 O 1 1
1 C 2 2
2 C 2 2
3 N 1 0
4 C 0 0
5 N 1 1
6 N 1 0
7 C 2 2
8 C 2 2
9 O 0 0
10 C 0 0
11 C 1 1
12 C 1 1
13 C 0 0
14 Cl 0 0
15 C 1 1
16 C 0 0
17 Cl 0 0
18 C 0 0
19 C 1 1
20 C 1 1
21 C 1 1
22 C 1 1
23 C 0 0

kekulized form (count_implicit_h(s2)):

0 O 1 1
1 C 2 2
2 C 2 2
3 N 0 0
4 C 0 0
5 N 1 1
6 N 0 0
7 C 2 2
8 C 2 2
9 O 0 0
10 C 0 0
11 C 1 1
12 C 1 1
13 C 0 0
14 Cl 0 0
15 C 1 1
16 C 0 0
17 Cl 0 0
18 C 0 0
19 C 1 1
20 C 1 1
21 C 1 1
22 C 1 1
23 C 0 0
$\endgroup$
1
  • 1
    $\begingroup$ Wow! Thanks! Its sad to know that pysmiles had a bug all along :( $\endgroup$
    – nir shahar
    Jul 31 at 11:28
7
$\begingroup$

Speculating a comparison of the outputs by both programs may help, I let RDKit process to SMILES to yield a representation including the atom indices. Derived from the cookbook, the MWE for the Jupyter notebook was

from rdkit import Chem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
from rdkit.Chem import rdDepictor
from rdkit.Chem.Draw import rdMolDraw2D

from IPython.display import SVG
IPythonConsole.ipython_useSVG=True  


def mol_with_atom_index(mol):
    for atom in mol.GetAtoms():
        atom.SetAtomMapNum(atom.GetIdx())
    return mol


mol = Chem.MolFromSmiles("OCCn2c(=N)n(CCOc1ccc(Cl)cc1Cl)c3ccccc23")
mol = mol_with_atom_index(mol)
mc = Chem.Mol(mol.ToBinary())


drawer = rdMolDraw2D.MolDraw2DSVG(450, 200) 
drawer.DrawMolecule(mc)
drawer.FinishDrawing()

svg = drawer.GetDrawingText()
display(SVG(svg.replace('svg:','')))

In a second step, I requested pysmiles to list for non-H atoms atom index, atom symbol (first two columns) and number of hydrogen bound (trailing column).

from pysmiles import read_smiles

smiles = "OCCn2c(=N)n(CCOc1ccc(Cl)cc1Cl)c3ccccc23"
molecule = read_smiles(smiles)

pysmiles_list = zip(molecule.nodes(data="element"), molecule.nodes(data="hcount"))
for element in pysmiles_list:
    print(element)

By visual inspection, RDKit's annotated formula and pysmiles' listing agree with each other except for the endocyclic nitrogens (N:3 and N:6): enter image description here

((0, 'O'), (0, 1))
((1, 'C'), (1, 2))
((2, 'C'), (2, 2))
((3, 'N'), (3, 1))
((4, 'C'), (4, 0))
((5, 'N'), (5, 1))
((6, 'N'), (6, 1))
((7, 'C'), (7, 2))
((8, 'C'), (8, 2))
((9, 'O'), (9, 0))
((10, 'C'), (10, 0))
((11, 'C'), (11, 1))
((12, 'C'), (12, 1))
((13, 'C'), (13, 0))
((14, 'Cl'), (14, 0))
((15, 'C'), (15, 1))
((16, 'C'), (16, 0))
((17, 'Cl'), (17, 0))
((18, 'C'), (18, 0))
((19, 'C'), (19, 1))
((20, 'C'), (20, 1))
((21, 'C'), (21, 1))
((22, 'C'), (22, 1))
((23, 'C'), (23, 0))

A similar listing of atom index, atom type, and number of hydrogens bond to the non-H may be achieved iterating with RDKit over the non-H atoms

for atom in mol.GetAtoms():
    print("{:2} {:2} {}".format(atom.GetIdx(), atom.GetSymbol(),
                                atom.GetTotalNumHs()))

to yield

 0 O  1
 1 C  2
 2 C  2
 3 N  0
 4 C  0
 5 N  1
 6 N  0
 7 C  2
 8 C  2
 9 O  0
10 C  0
11 C  1
12 C  1
13 C  0
14 Cl 0
15 C  1
16 C  0
17 Cl 0
18 C  0
19 C  1
20 C  1
21 C  1
22 C  1
23 C  0

As outlined in the answer by @rapelpy, there is a problem with pysmiles which the maintainer(s) might rectify in future. The computations in this answer were performed with pysmiles (version 1.0.1, Jun 25, 2020), and RDKit as packaged for Debian 11 (testing) (version 2020.09.4).

$\endgroup$
6
  • $\begingroup$ Hi, thanks for the help! I didn't quite get where I see the number of hydrogens that RDkit calculated here. Where is it shown? $\endgroup$
    – nir shahar
    Jul 29 at 20:21
  • $\begingroup$ Its only implicitly present, in the structure drawing. I wasn't sure if the atom indices used by RDKit and PySmiles would align; the map annotated by atom indices side by side to the listing now suggest they do. $\endgroup$
    – Buttonwood
    Jul 29 at 20:51
  • $\begingroup$ Hmm. How can I be sure that they align for all my data, without checking for a few examples and just hoping it won't do problems? $\endgroup$
    – nir shahar
    Jul 29 at 20:56
  • $\begingroup$ The Phonopy room has now been changed into the "python-based software room" since it wasn't very active before, and questions like this which are about a python software like PySmiles not matching with the output of RDKit, can benefit from having a room when the discussions in comments get long. $\endgroup$ Jul 29 at 21:56
  • 2
    $\begingroup$ @nirshahar continuation at chat.stackexchange.com/transcript/message/58750154#58750154 $\endgroup$
    – Buttonwood
    Jul 29 at 22:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.