6. Reactor

1. 2019, 2020 Dr. Ramil Nugmanov;
1. 2019 Dr. Timur Madzhidov; Ravil Mukhametgaleev
1. 2020-2022 Adelia Fatykhova
1. 2022 Valentina Afonina

Installation instructions of CGRtools package information and tutorial’s files see on https://github.com/cimm-kzn/CGRtools

NOTE: Tutorial should be performed sequentially from the start. Random cell running will lead to unexpected results.

[1]:

import pkg_resources
if pkg_resources.get_distribution('CGRtools').version.split('.')[:2] != ['4', '1']:
    print('WARNING. Tutorial was tested on 4.1 version of CGRtools')
else:
    print('Welcome!')

Welcome!

[2]:

# load data for tutorial
from pickle import load
from traceback import format_exc

with open('reactions.dat', 'rb') as f:
    reactions = load(f) # list of ReactionContainer objects

r1 = reactions[0] # reaction
m6 = r1.reactants[1]
m6copy = m6.copy()
m6copy.atom(5)._Core__isotope = 13

Reactor objects stores single transformation and can apply it to molecules or CGRs.

Transformations is ReactionContainer object which in reactant side consist of query for matching group and in product side patch for updating matched atoms and bonds

[3]:

from CGRtools import CGRReactor, Reactor   # import of Reactor
from CGRtools.containers import * # import of required objects
from CGRtools.containers.bonds import DynamicBond

6.1. Products generation

Reactor works similar to ChemAxon Reactions enumeration.

Example here presents application of it to create esters from acids.

First we need to construct carboxy group matcher query. Then, ether group need to be specified.

Atom numbers in query and patch should be mapped to each other. The same atoms should have same numbers.

[4]:

acid = QueryContainer()                       # this query matches acids. Use construction possibilities.
acid.add_atom('C', neighbors=3)   # add carboxyl carbon. Hybridization is irrelevant here
acid.add_atom('O', neighbors=1)   # add hydroxyl oxygen. Hybridization is irrelevant here
acid.add_atom('O')                                # add carbonyl oxygen. Number of neighbors is irrelevant here.
acid.add_bond(1, 2, 1) # create single bond between carbon and hydroxyl oxygen
acid.add_bond(1, 3, 2) # create double bond
print(acid)
acid.clean2d()
acid

[O]=[C;D3]-[O;D1]

[4]:

[5]:

methyl_ester = QueryContainer()  # create patch - how carboxyl group should be changed. We write methylated group
methyl_ester.add_atom('C', 1) # second argument is predefined atom mapping. Notice that mapping corresponds...
methyl_ester.add_atom('O', 2) # ... to order in already created acid group. Atom 2 is released water.
methyl_ester.add_atom('O', 4)
methyl_ester.add_atom('O', 3)
methyl_ester.add_atom('C', 5)
methyl_ester.add_bond(1, 4, 1)
methyl_ester.add_bond(1, 3, 2)
methyl_ester.add_bond(4, 5, 1)
# No bond between atom 1 and atom 2. This bond will be broken.
methyl_ester.clean2d()
methyl_ester

[5]:

[6]:

m6 # acid

[6]:

[7]:

template = ReactionContainer([acid], [methyl_ester]) # merge query and patch in template, which is ReactionContainer
reactor = CGRReactor(template)                        # Reactor is initialized
reacted_acid = next(reactor(m6))                            # application of Reactor to molecule

[8]:

reacted_acid.clean2d() # calculate coordinates
reacted_acid       # desired methylated ester have been generated

[8]:

One can notice presence of separate oxygen (water) and ester group.

The second group can substituted by calling reactor on observed product.

[9]:

second_stage = next(reactor(reacted_acid)) # apply transformation on product of previous transformation
second_stage.clean2d() #  recalculate coordinates for correct drawing
second_stage

[9]:

second_stage has 3 components in a single MoleculeContainer object. We can split it into individual molecules and place all molecules into ReactionContainer object. Since in CGRtools atom-to-atom mapping corresponds to numbering of atoms in molecules, the resulting product has AAM according to the rule applied. Thus, reaction has correct AAM and nothing special should be made to keep or find it.

[10]:

products = second_stage.split() # split product into individual molecules
react = ReactionContainer([m6], products) # unite reagent and product into reaction.
react

[10]:

For multicomponent reactions one can merge molecules of reactants into single MoleculeContainer object and apply reactor on it.

It is possible to generate all available products in case that molecule has several groups matching the query.

[11]:

m6copy

[11]:

[12]:

enums = set()              # the set enums is used to select structurally diverse products
for m in reactor(m6copy): # limit=0 is enumeration of all possible products by reactor
    print(m)                       # print signatures for observed molecules. Notice presence of water as component of product
    m.clean2d()         # recalculate coordinates
    enums.update(m.split())        # split product into separate molecules
enums = list(enums)                # set of all resulting molecules

O.O=C([13C](O)=O)OC
O=C(C(OC)=O)O.O

Let’s have a look at molecules in set. Note to lost isotope mark.

[13]:

enums[0]

[13]:

[14]:

enums[1]

[14]:

[15]:

enums[2]

[15]:

6.2. MetaReactions (reactions on CGRs).

Reactor could be applied to CGR to introduce changes into reaction.

6.2.1. Example of atom-to-atom mapping fixing.

[16]:

reactions[1] # reaction under study

[16]:

[17]:

cgr = ~reactions[1] # generate reaction CGR
print(cgr)
cgr.clean2d()
cgr

c1cc(ccc1[N+]([->=][O->0])[=>-][O0>-])O[.>-]C(C)[->.]I

[17]:

[18]:

cgr.centers_list # reaction has two reaction centers. [10,11,12] - pseudo reaction appeared due to AAM error

[18]:

((2, 3, 13), (10, 11, 12))

Reaction has AAM error in nitro-group

Lets try to use Reactor for AAM fixing

[19]:

nitro = QueryCGRContainer() # construct query for invalid reaction center - CGR of wrongly mapped nitro-group
nitro.add_atom('N', charge=1, p_charge=1) # atom 1
nitro.add_atom('O', charge=0, p_charge=-1) # atom 2. Notice that due to AAM error charge was changed
nitro.add_atom('O', charge=-1, p_charge=0) # atom 3. Notice that due to AAM error charge was changed
nitro.add_atom('C') # atom 4

nitro.add_bond(1, 2, DynamicBond(2, 1)) # bond between atoms 1 and 2. Due to AAM error bond is dynamic ('2>1' type)
nitro.add_bond(1, 3, DynamicBond(1, 2)) # bond between atoms 1 and 3. Due to AAM error bond is dynamic ('1>2' type)
nitro.add_bond(1, 4, 1) # ordinary bond
print(nitro)
# this query matches reaction center in CGR appeared due to AAM error.
nitro.clean2d()
nitro

[C][N+]([=>-][O0>-])[->=][O->0]

[19]:

[20]:

nitro < cgr # query matches CGR of reaction with error.

[20]:

True

[21]:

valid_nitro = QueryCGRContainer() # construct nitro group without dynamic atoms. Notice that atom order should correspond object nitro
valid_nitro.add_atom('N', charge=1, p_charge=1) # ordinary N atom
valid_nitro.add_atom('O', charge=-1, p_charge=-1) # ordinary negatively charged oxygen atom
valid_nitro.add_atom('O')                            # ordinary oxygen atom
valid_nitro.add_atom('C') # atom 4

valid_nitro.add_bond(1, 2, 1) # ordinary single bond
valid_nitro.add_bond(1, 3, 2) # ordinary double bond
valid_nitro.add_bond(1, 4, 1) # ordinary bond
print(valid_nitro)
valid_nitro.clean2d()
valid_nitro

[O]=[N+]([C])[O-]

[21]:

Now time to prepare and apply Template to CGR based on reaction with incorrect AAM.

Template is Reaction container with query in reactants and patch in products

[22]:

template = ReactionContainer([nitro], [valid_nitro]) # template shows how wrong part of CGR is transformed into correct one.
print(template) # notice complex structure of query: CGR signature is given in braces, then >> and molecule signature
template

[C][N+]([=>-][O0>-])[->=][O->0]>>[O]=[N+]([C])[O-]

[22]:

Reactor class accept single template. Existence of dynamic bond in it is not a problem.

[23]:

reactor = CGRReactor(template)

Reactor object is callable and accept as argument molecule or CGR.

NOTE: fixed is new CGR object

[24]:

fixed = next(reactor(cgr)) # fix CGR

CGRreactor returns None if template could not be applied, otherwise patched structure is returned.

[25]:

print(fixed)
fixed

c1cc(ccc1[N+](=O)[O-])O[.>-]C(C)[->.]I

[25]:

One can see that nitro group has no dynamic bonds any more. CGR corresponds only to substitution.

[26]:

fixed.centers_list # reaction center appeared due to AAM error before does not exist. Only 1 reaction center is found

[26]:

((2, 3, 13),)

6.2.2 Reaction transformation

Example of E2 to SN2 transformation.

E2 and SN2 are concurrent reactions. We can easily change reaction center of E2 reaction to SN2. It could be achieved by substitution of reaction center corresponding to double bond formation in E2 reaction by the one corresponding to formation of new single bond with base as in SN2.

[27]:

from CGRtools.files import MRVRead
from io import StringIO

[28]:

e2 = next(MRVRead('e2.mrv')) # read E2 reaction from ChemAxon MRV file
e2

[28]:

[29]:

# create CGR query for E2 reaction side
e2query = QueryCGRContainer()
e2query.add_atom('C', 1) # create carbon with mapping number 1
e2query.add_atom('C', 2) # create carbon with mapping number 2
# addition of iodine atom
e2query.add_atom('I', 3, neighbors=1, p_neighbors=0, charge=0, p_charge=-1)
# addition of OH- or RO- groups
e2query.add_atom('O', 4, neighbors=[0, 1], p_neighbors=[0, 1], charge=-1, p_charge=0)

e2query.add_bond(1, 2, DynamicBond(1, 2)) # bond between two carbon corresponds to formation of double from single
e2query.add_bond(1, 3, DynamicBond(1)) # bond between carbon and halogen breaks in E2 reaction
print(e2query) # it is CGR of E2 reaction center
e2query.clean2d()
e2query

[C][->=][C][->.][I;1>0;0>-].[O;01>01;->0]

[29]:

[30]:

e2_cgr = ~e2 # compose reaction into CGR
e2_cgr

[30]:

[31]:

e2query < e2_cgr # E2 CGR pattern works!

[31]:

True

[32]:

# create patch creating SN2 reaction. Notice that ordering of atoms correspond to that of E2 CGR query
sn2patch = QueryCGRContainer()
sn2patch.add_atom('C', 1) # save atom unchanged.
sn2patch.add_atom('C', 2) # it is central atom.
sn2patch.add_atom('I', 3, charge=0, p_charge=-1)
sn2patch.add_atom('O', 4, charge=-1, p_charge=0)

sn2patch.add_bond(1, 2, 1) # set carbon - carbon single bond that is unchanged in SN2 reaction
sn2patch.add_bond(1, 3, DynamicBond(1)) # this bond is broken in SN2 reaction
sn2patch.add_bond(1, 4, DynamicBond(None, 1)) # it corresponds to formation of bond O(S)-C bond in SN2 reaction
print(sn2patch)
sn2patch.clean2d()
sn2patch

[I0>-][->.][C]([.>-][O->0])[C]

[32]:

[33]:

reactor = CGRReactor(ReactionContainer([e2query], [sn2patch])) # create template and pass it to Reactor
sn2_cgr = next(reactor(e2_cgr)) # apply Reactor on E2 reaction

[34]:

print(sn2_cgr)
sn2_cgr

C(CC)(C)([->.][I0>-])[.>-][O->0].[Na+]

[34]:

[35]:

# decompose CGR into reaction
sn2 = ReactionContainer.from_cgr(sn2_cgr)
sn2.clean2d()
sn2

[35]: