Different behaviour of molecules returned by RGroupDecompose between RDKit v2022.9.5 & 2023.3

[LindeSchoenmaker]

I am using rdkit to transfer R-groups from one molecule to another.

My workflow is as follows;

1. decompose a pair of molecules into a common core and R-groups
2. pick a combination of R-groups
3. attach those to the common core

I noticed that the molecules returned by RGroupDecompose (I'm using asSmiles = False), behave differently between the different versions.

In version 2022.9.5 I'm able to add two R-groups to a positively charged primary amine (C[NH2+][*:1]). For version 2023.3 the same code leads to a valence error.

I noticed that the cores returned for v2022.9.5 are all identical, whereas for v2023.3 the number of explicit Hs in the SMILES is the same as for the parent molecule.

Below I put the code to reproduce my results.

I'm currently post-processing my molecules to fix the number of hydrogens by counting the number of bonds and comparing that to the allowed valence. However, it would be preferable if the change in returned molecules (from v.2023.3 onwards) could be prevented or the molecules could be transformed 'back' to what they were like in the earlier versions.

import rdkit

print(rdkit.__version__)
import pandas as pd
from rdkit import Chem
from rdkit.Chem import rdRGroupDecomposition as rdRGD
from rdkit.Chem.rdchem import EditableMol

smiles = ['[NH3+]C[C@H]1CC[C@@H]2[C@H](O1)C1=CC(C(C)(C)C)=CC=C1N[C@H]2C1=CC=CC=C1', 'CC[N@H+](CCO)C[C@H]1CC[C@@H]2[C@H](O1)C1=CC(C(C)(C)C)=CC=C1N[C@H]2C1=CC=CC=C1']
mols = [Chem.MolFromSmiles(x) for x in smiles]

res = Chem.rdFMCS.FindMCS(mols,
                        matchValences=True,
                        ringMatchesRingOnly=True,
                        completeRingsOnly=True,
                        timeout=2)
core = Chem.MolFromSmarts(res.smartsString)

res, _ = rdRGD.RGroupDecompose([core],
                                mols,
                                asSmiles=False,
                                asRows=False)
df_rgroup = pd.DataFrame(res)
mol_to_weld = Chem.CombineMols(df_rgroup['Core'][0], df_rgroup['R1'][1])

# transfer atom map numbers and remove dummy atoms
idxs = []
for atom in mol_to_weld.GetAtoms():
    map_num = atom.GetAtomMapNum()
    if map_num > 0:
        idxs.append((atom.GetAtomMapNum(), atom))

for idx, atom in idxs:
    nbr_1 = [x.GetOtherAtom(atom) for x in atom.GetBonds()][0]
    nbr_1.SetAtomMapNum(idx)

new_mol = Chem.DeleteSubstructs(mol_to_weld, Chem.MolFromSmarts('[#0]'))

# get index of atoms to connect

idxs = []
for atom in new_mol.GetAtoms():
    map_num = atom.GetAtomMapNum()
    if map_num > 0:
        idxs.append(atom.GetIdx())

# add bonds
em = EditableMol(new_mol)
em.AddBond(idxs[0], idxs[1], order=Chem.BondType.SINGLE)
em.AddBond(idxs[0], idxs[2], order=Chem.BondType.SINGLE)
mol_new = em.GetMol()

[bp-kelley]

I can't speak to the current change in behavior, however Chem.molzip was designed to recombine the results of a decomposition. There is one caveat that the Core must come first otherwise molzip gets confused.

The one other change is that when using rdFMCS, the decomposition behaves better when using res.queryMol as opposed to res.smartsString. This simplifies your example to the following (similar to what we use in the FreeWilson contribution)

smiles = ['[NH3+]C[C@H]1CC[C@@H]2[C@H](O1)C1=CC(C(C)(C)C)=CC=C1N[C@H]2C1=CC=CC=C1', 'CC[N@H+](CCO)C[C@H]1CC[C@@H]2[C@H](O1)C1=CC(C(C)(C)C)=CC=C1N[C@H]2C1=CC=CC=C1']
mols = [Chem.MolFromSmiles(x) for x in smiles]

res = rdFMCS.FindMCS(mols,
                        matchValences=True,
                        ringMatchesRingOnly=True,
                        completeRingsOnly=True,
                        timeout=2)

res, _ = rdRGD.RGroupDecompose([res.queryMol],
                                mols,
                                asSmiles=False,
                                asRows=False)
# for molzip, core must come first when combining RGroupDecompose objects
#. normally does by default
assert list(res.keys())[0] == 'Core'

import itertools
# run through all combinations of cores and rgroups
for i, rgroups in enumerate(itertools.product(*res.values())):
    # make a molecule with Core+R1+R2+R3...
    mol = rgroups[0]
    for rg in rgroups[1:]:
        mol = Chem.CombineMols(mol, rg)
    # Zip the molecule together
    res = Chem.molzip(mol)
    print(i, Chem.MolToSmiles(res))

[bp-kelley]

There is also a benefit to molzip in that it handles stereochemistry as well.

[LindeSchoenmaker]

Thank you for putting Chem.molzip on my radar. This indeed does not exceed the allowed valence.