How to change (or remap) atom index (ordering) of a molecule after canonicalization?

[qcxia20]

Hi all,
I have a question on that if there are any methods to map the atom indexes (I mean, ordering in the molecule) of the molecule before and after canonicalization.
I'm doing the canonicalization through the way below. I found the atom indexes (ordering) changed because that the canonicalization step has changed the raw orientation of mol and then the newly generated mol canonical_mol is assigned with different atom ordering.

mol = Chem.MolFromMol2File(mol2file)
canonical_mol = Chem.MolFromSmiles(Chem.MolToSmiles(mol))

More specifically, my problem is that I have to locate the same groups of torsional bonds on the both two molecules with one-by-one mapping. The location of torsional bonds are determined through the way below. And because 1) the orientation and atom ordering are different in mol and canonical_mol and 2) I cannot get the specific torsion SMARTS through this way, the list all_rot which contains all rotatable bonds in the molecule cannot be mapped one-by-one directly. (Say, rotatable bond 1 in mol can actually be rotatable bond 4 in canonical_mol and vice versa.)

from rdkit.Chem.Lipinski import RotatableBondSmarts
all_rot = mol.GetSubstructMatches(RotatableBondSmarts)

So any specific ideas on solving my specific problem are also welcome and appreciated.

Thanks!

[ptosco]

@qcxia20 MolFromSmiles adds two private properties to the molecule it is called on, _smilesAtomOutputOrder and _smilesBondOutputOrder, which contain the original atom index before canonicalization. See below:

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

IPythonConsole.drawOptions.addAtomIndices = True

mol = Chem.MolFromSmiles("OC(COC1=C2C=CC=CC2=CC=C1)CNC(C)C")
mol

canonical_mol = Chem.MolFromSmiles(Chem.MolToSmiles(mol))
canonical_mol

list(mol.GetPropNames(includePrivate=True, includeComputed=True))
['__computedProps',
 'numArom',
 '_StereochemDone',
 '_smilesAtomOutputOrder',
 '_smilesBondOutputOrder']
list(map(int, mol.GetProp("_smilesAtomOutputOrder")[1:-2].split(",")))
[17, 16, 18, 15, 14, 1, 0, 2, 3, 4, 13, 12, 11, 10, 9, 8, 7, 6, 5]

So atom 0 was atom 17 before canonicalization, 1 was 16, 2 was 18, etc.

[qcxia20]

Thanks Paolo! This approach works.

[tsn2027]

can I have a follow-up question about this? It seems like _smilesAtomOutputOrder does not show up in rdkit's version 2023.03.2 anymore. May I ask how I can remap atom index after canonicalization with a new version? Thanks.

[greglandrum]

@tsn2027 I don't have a local environment with an install of 2023.03.2, but we have definitely not removed that functionality from the RDKit:

In [9]: m = Chem.MolFromSmiles('OCC')

In [10]: Chem.MolToSmiles(m)
Out[10]: 'CCO'

In [11]: m.GetProp('_smilesAtomOutputOrder')
Out[11]: '[2,1,0,]'

In [12]: rdkit.__version__
Out[12]: '2023.09.5'

Could it be that you forgot to call Chem.MolToSmiles() on the molecule first?