atom mapping in reaction

[manajit-das]

I am creating reaction smiles using AllChem.ReactionFromSmarts as shown here: https://gist.github.com/greglandrum/8028f785fe3ab0162991b22f1f4a2653. I want the mapping of atom index between reactant and product. The following function https://github.com/rdkit/rdkit/issues/1269#issuecomment-278205412

def make_rxns(source_rxn, reactants):
    product_sets = source_rxn.RunReactants(reactants)
    new_rxns = []
    for pset in product_sets:
        new_rxn = AllChem.ChemicalReaction()
        for react in reactants:
            react = Chem.Mol(react)
            for a in react.GetAtoms():
                a.SetIntProp('molAtomMapNumber', a.GetIdx())
            new_rxn.AddReactantTemplate(react)
        for prod in pset:
            for a in prod.GetAtoms():
                a.SetIntProp('molAtomMapNumber', int(a.GetProp('react_atom_idx')))
            new_rxn.AddProductTemplate(prod)
        new_rxns.append(new_rxn)
    return new_rxns

gives a reasonable answer but when I try a reaction with two reactants, every reactant's index start from 0, hence there are overlapping index between reactants. For example:

mol1 = Chem.MolFromSmiles('Clc1ccccc1')
mol2 = Chem.MolFromSmiles('C[P+]([Pd])(C)C')
sr = AllChem.ReactionFromSmarts('[c:3][Cl:1].[Pd:2]>>[Pd:2]([Cl:1])[c:3]')

for r in make_rxns(sr, [mol1, mol2]):
    print(AllChem.ReactionToSmiles(r))
    rxn = AllChem.ReactionToSmiles(r)

[CH3:0][P+:1]([Pd:2])([CH3:3])[CH3:4].[Cl:0][c:1]1[cH:2][cH:3][cH:4][cH:5][cH:6]1>>[CH3:0][P+:1]([Pd:2]([Cl:0])[c:1]1[cH:2][cH:3][cH:4][cH:5][cH:6]1)([CH3:3])[CH3:4]

So, the indexes 0, 1, 2, 3, 4 are repeated in the reactants. I would like to have unique ids for every atom which are mapped to the products.
Thank you.

[greglandrum]

It's tricky (possible, but tricky) to do this correctly. But there's an easy workaround.

I think your best bet here is to combine the reactants into a single molecule and add parentheses to the reaction SMARTS to group the two reactants into one:

mol1 = Chem.MolFromSmiles('Clc1ccccc1')
mol2 = Chem.MolFromSmiles('C[P+]([Pd])(C)C')
rmol = Chem.CombineMols(mol1,mol2)
sr = AllChem.ReactionFromSmarts('([c:3][Cl:1].[Pd:2])>>[Pd:2]([Cl:1])[c:3]')
for r in make_rxns(sr, [rmol]):
    smi = AllChem.ReactionToSmiles(r)
    print(smi)

This outputs:

([CH3:7][P+:8]([Pd:9])([CH3:10])[CH3:11].[Cl:0][c:1]1[cH:2][cH:3][cH:4][cH:5][cH:6]1)>>[Cl:0][Pd:9]([c:1]1[cH:2][cH:3][cH:4][cH:5][cH:6]1)[P+:8]([CH3:7])([CH3:10])[CH3:11]

that is almost what you want, but it still has the parentheses around the reactants. You can remove those with a regular expression:

for r in make_rxns(sr, [rmol]):
    smi = AllChem.ReactionToSmiles(r)
    smi =  re.sub(r'^\((.*)\)>',r'\1>',smi)
    print(smi)

This outputs:

[CH3:7][P+:8]([Pd:9])([CH3:10])[CH3:11].[Cl:0][c:1]1[cH:2][cH:3][cH:4][cH:5][cH:6]1>>[Cl:0][Pd:9]([c:1]1[cH:2][cH:3][cH:4][cH:5][cH:6]1)[P+:8]([CH3:7])([CH3:10])[CH3:11]

I think this ought to work

[manajit-das]

Thank you so much:) it works.

[greglandrum]

awesome!