Merge pull request #47 from ykevu/atom-map-reactants

Allow for custom atom mapping in reactants when running reaction from SMARTS

rdchiral/initialization.py

@@ -93,12 +93,13 @@ class rdchiralReactants(object):
     Args:
         reactant_smiles (str): Reactant SMILES string
     '''
-    def __init__(self, reactant_smiles):
+    def __init__(self, reactant_smiles, custom_reactant_mapping=False):
         # Keep original smiles, useful for reporting
         self.reactant_smiles = reactant_smiles
+        self.custom_mapping = custom_reactant_mapping
 
         # Initialize into RDKit mol
-        self.reactants = initialize_reactants_from_smiles(reactant_smiles)
+        self.reactants = initialize_reactants_from_smiles(reactant_smiles, custom_reactant_mapping)
 
         # Set mapnum->atom dictionary
         # all reactant atoms must be mapped after initialization, so this is safe
@@ -107,7 +108,7 @@ def __init__(self, reactant_smiles):
 
         # Create copy of molecule without chiral information, used with
         # RDKit's naive runReactants
-        self.reactants_achiral = initialize_reactants_from_smiles(reactant_smiles)
+        self.reactants_achiral = initialize_reactants_from_smiles(reactant_smiles, custom_reactant_mapping)
         [a.SetChiralTag(ChiralType.CHI_UNSPECIFIED) for a in self.reactants_achiral.GetAtoms()]
         [(b.SetStereo(BondStereo.STEREONONE), b.SetBondDir(BondDir.NONE)) \
             for b in self.reactants_achiral.GetBonds()]
@@ -167,7 +168,7 @@ def initialize_rxn_from_smarts(reaction_smarts):
 
     return rxn
 
-def initialize_reactants_from_smiles(reactant_smiles):
+def initialize_reactants_from_smiles(reactant_smiles, custom_reactant_mapping):
     '''Initialize RDKit molecule from SMILES string
 
     Args:
@@ -183,7 +184,8 @@ def initialize_reactants_from_smiles(reactant_smiles):
     # To have the product atoms match reactant atoms, we
     # need to populate the map number field, since this field
     # gets copied over during the reaction via reactant_atom_idx.
-    [a.SetAtomMapNum(i+1) for (i, a) in enumerate(reactants.GetAtoms())]
+    if not custom_reactant_mapping:
+        [a.SetAtomMapNum(i+1) for (i, a) in enumerate(reactants.GetAtoms())]
     if PLEVEL >= 2: print('Initialized reactants, assigned map numbers, stereochem, flagpossiblestereocenters')
     return reactants
 

rdchiral/main.py

@@ -76,7 +76,7 @@
 
 '''
 
-def rdchiralRunText(reaction_smarts, reactant_smiles, **kwargs):
+def rdchiralRunText(reaction_smarts, reactant_smiles, custom_reactant_mapping=False, **kwargs):
     '''Run from SMARTS string and SMILES string. This is NOT recommended
     for library application, since initialization is pretty slow. You should
     separately initialize the template and molecules and call run()
@@ -90,7 +90,7 @@ def rdchiralRunText(reaction_smarts, reactant_smiles, **kwargs):
         list: List of outcomes from `rdchiralRun`
     '''
     rxn = rdchiralReaction(reaction_smarts)
-    reactants = rdchiralReactants(reactant_smiles)
+    reactants = rdchiralReactants(reactant_smiles, custom_reactant_mapping)
     return rdchiralRun(rxn, reactants, **kwargs)
 
 def rdchiralRun(rxn, reactants, keep_mapnums=False, combine_enantiomers=True, return_mapped=False):

test/test_atom_mapping.py

@@ -0,0 +1,42 @@
+import os, sys, json
+sys.path = [os.path.dirname(os.path.dirname((__file__)))] + sys.path 
+
+from rdchiral.main import rdchiralReaction, rdchiralReactants, rdchiralRunText, rdchiralRun
+from rdkit import Chem
+
+with open(os.path.join(os.path.dirname(__file__), 'test_atom_mapping_cases.json'), 'r') as fid:
+    test_cases = json.load(fid)
+
+all_passed = True
+
+def canonicalize_outcomes(outcomes):
+    ''' Convert all SMILES in a list of outcomes to the canonical form '''
+    return list(map(lambda x: Chem.CanonSmiles(x), outcomes))
+
+for i, test_case in enumerate(test_cases):
+
+    print('\n# Test {:2d}/{}'.format(i+1, len(test_cases)))
+
+    # Initialize with test case SMILES/SMARTS
+    reaction_smarts = test_case['smarts']
+    reactant_smiles = test_case['smiles']
+    reactants = rdchiralReactants(reactant_smiles, custom_reactant_mapping=True)
+    expected = canonicalize_outcomes(test_case['expected'])
+
+    # Test rdchiralRunText
+    if canonicalize_outcomes(rdchiralRunText(reaction_smarts, reactant_smiles, custom_reactant_mapping=True, keep_mapnums=True)) == expected:
+        print('    from text: passed')
+    else:
+        print('    from text: failed')
+        all_passed = False
+
+    # Pre-initialize & repeat with rdChiralRun
+    rxn = rdchiralReaction(reaction_smarts)
+    if all(canonicalize_outcomes(rdchiralRun(rxn, reactants, keep_mapnums=True)) == expected for j in range(3)):
+        print('    from init: passed')
+    else:
+        print('    from init: failed')
+        all_passed = False
+
+all_passed = 'All passed!' if all_passed else 'Failed!'
+print('\n# Final result: {}'.format(all_passed))

test/test_atom_mapping_cases.json

@@ -0,0 +1,50 @@
+[
+    {
+        "smarts": "[C;R:10][O;R:2][Xe;R:3][O;R:4][C;R:5]>>[C:10][O:2].[O;R:4][C;R:5]",
+        "smiles": "[C:10]1[C:9][O:8][Xe:7][O:6]1",
+        "expected": ["[C:10]([C:9][OH:8])[OH:6]"],
+        "description": "Test custom atom mapping in reactant"
+    },
+	{
+		"smarts": "[C:1](=[O:3])[OH:2]>>[C:1](=[O:3])[O:2]CC",
+		"smiles": "[OH:10][C:11](=[O:12])[CH2:13][CH2:14][CH2:15][CH2:16][CH2:17][CH3:18]",
+		"expected": ["[CH3:18][CH2:17][CH2:16][CH2:15][CH2:14][CH2:13][C:11](=[O:12])[O:10][CH2:900][CH3:901]"],
+		"description": "Testing achiral transformations with chiral molecules (not in template): Preparing a carboxylic acid from hydrolysis of an ethyl ester"
+	},
+	{
+		"smarts": "[C:1](=[O:3])[OH:2]>>[C:1](=[O:3])[O:2]CC",
+		"smiles": "[OH:10][C:11](=[O:12])[CH2:13][CH2:14][CH2:15][CH2:16][C@H:17]([Cl:18])[C:19]",
+		"expected": ["[CH3:901][CH2:900][O:10][C:11](=[O:12])[CH2:13][CH2:14][CH2:15][CH2:16][C@@H:17]([C:19])[Cl:18]"],
+		"description": "Testing achiral transformations with chiral molecules (not in template): Preparing a carboxylic acid from hydrolysis of an ethyl ester"
+	},
+	{
+		"smarts": "[C:4][C:1](=[O:3])[OH:2]>>[C:4][C:1](=[O:3])[O:2]CC",
+		"smiles": "[OH:10][C:11](=[O:12])[C@H:13]([Cl:14])[CH3:15]",
+		"expected": ["[CH3:901][CH2:900][O:10][C:11](=[O:12])[C@@H:13]([CH3:15])[Cl:14]"],
+		"description": "Testing achiral transformations with chiral molecules (partially in template, but auxiliary): Preparing a carboxylic acid from hydrolysis of an ethyl ester"
+	},
+	{
+		"smarts": "[C:1][CH:2]([CH3:3])[O:4][C:5]>>[C:1][CH:2]([CH3:3])[OH:4].O[C:5]",
+		"smiles": "[CH3:10][CH2:11][CH2:12][CH2:13][C@@H:14]([O:15][CH2:16][CH2:17])[CH3:18]",
+		"expected": [],
+		"description": "Testing achiral transformations with chiral molecules (fully in template): Alkylation reaction with unspecified chirality, template could have specified"
+	},
+	{
+		"smarts": "[C:1][C@H:2]([CH3:3])[I:4]>>[C:1][C@@H:2]([CH3:3])Br",
+		"smiles": "[CH3:10][CH2:11][CH2:12][CH2:13][CH:14]([I:15])[CH3:16]",
+		"expected": [],
+		"description": "Testing chiral transformations with achiral molecules: SN2 with inversion of a tetrahedral center"
+	},
+	{
+		"smarts": "[C:1]/[CH:2]=[CH:3]\\[C:4]>>[C:1][CH0:2]#[CH0:3][C:4]",
+		"smiles": "[CH3:10][CH2:11][CH2:12]/[CH:13]=[CH:14]\\[CH2:15][CH3:16]",
+		"expected": ["[CH3:16][CH2:15][C:14]#[C:13][CH2:12][CH2:11][CH3:10]"],
+		"description": "Testing chiral transformations with chiral molecules: Reaction expects cis double bond: Molecule has explicit cis double bond: explicit cis"
+	},
+	{
+		"smarts": "[C:1](=[O:3])[O:2][C:4]>>[C:1](=[O:3])[OH:2].O[C:4]",
+		"smiles": "[CH2:10]1[C:11](=[O:12])[O:13][CH2:14][CH2:15][CH2:16]1",
+		"expected": ["[O:12]=[C:11]([OH:13])[CH2:10][CH2:16][CH2:15][CH2:14][OH:900]"],
+		"description": "Accidental fragmentation: intramolecular esterification"
+	}
+]