Unexpected HasReactionSubstructMatch behavior?

[thorlucas]

Let's say I'm defining a reaction where any R-CC group can be converted to R-CN group. As I understand it, I would write this in reaction SMARTS like this: CC>>CN.

Then, let's say we're querying for some reaction that would produce the product CCCN. So in SMARTS we write >>CCCN.

Then we query:

rxn = Reactions.ReactionFromSmarts('CC>>CN')
q = Reactions.ReactionFromSmarts('>>CCCN')
Reactions.HasReactionSubstructMatch(rxn, q)
>>> False

Why is this false?

[greglandrum]

Hi @thorlucas,
HasReactionSubstructMatch(r1,r2) looks to see if r2 is a substruct of r1, not the other way around.
So you get a match for something like this:

In [6]: rxn = Reactions.ReactionFromSmarts('CC>>CN')
   ...: q = Reactions.ReactionFromSmarts('CCCC>>CCCN')
   ...: Reactions.HasReactionSubstructMatch(q,rxn)
Out[6]: True

In the example you provide I think what you're looking for is to determine whether or not a particular molecule could be a product of a reaction. You can do that with IsMoleculeProduct():

In [14]: rxn = Reactions.ReactionFromSmarts('CC>>CN')

In [15]: rxn.Initialize()
[05:41:24] reactant 0 has no mapped atoms.
[05:41:24] product 0 has no mapped atoms.

In [16]: rxn.IsMoleculeProduct(Chem.MolFromSmiles('CCCN'))
Out[16]: True

In [17]: rxn.IsMoleculeProduct(Chem.MolFromSmiles('CCCO'))
Out[17]: False

[greglandrum]

I'm not sure I follow. If we have the reaction 'CC>>CN' (a reaction which takes an alkane and converts it to an amine), then the reaction '>>CCN' (a reaction that produces an amine) should match it. Yet:

rxn = Reactions.ReactionFromSmarts('CC>>CN')
q = Reactions.ReactionFromSmarts('>>CCCN')

Reactions.HasReactionSubstructMatch(q, rxn)
>>> False

Reactions.HasReactionSubstructMatch(rxn, q)
>>> False

It looks like you are expecting the RDKit reaction code to do something other than what it does.
The function rdChemReactions.HasReactionSubstructMatch(r1,r2) looks to see if r2 is a substructure of r1.
From the RDKit's perspective that means that each reactant of r2 is a substructure of a reactant of r1 and each product of r2 is a substructure of a product of r1.

Additionally, it seems that this method does not respect atom mapping:

In [4]: def cmp(a, b):
   ...:     rxn = Reactions.ReactionFromSmarts(a)
   ...:     q = Reactions.ReactionFromSmarts(b)
   ...:     return Reactions.HasReactionSubstructMatch(q, rxn)

In [7]: cmp('[C:1][C:2]>>[C:1]N', '[C:1][C:2][C:3]>>[C:1][C:2]N')
Out[7]: True
# Makes sense

In [8]: cmp('[C:1][C:2]>>[C:1]N', '[C:1][C:2][C:3]>>[C:1][C:5]N')
Out[8]: True
# Should be false

In this reaction, we can see that the end product should have the amine attached to the original carbon C1. Yet, we can see that in the second example, I have it attached a totally different carbon, and it still hits.

That is correct. The substructure matching code currently ignores atom mapping information. This is something we should look into changing, particularly when doing substructure matches on reactions.

[greglandrum]

On the question about R groups: I'm not sure what you're trying to do. If you want to define a transformation which brominates a primary amine, you can do the following:

In [10]: rxn = rdChemReactions.ReactionFromSmarts('[CD2:1]-[OH]>>[C:1]-Br')

Notice that I've mapped the atoms so that I can actually use the reaction to produce products.

Here's a demonstration that this reaction actually does what I expect it to do:

In [11]: ps = rxn.RunReactants((Chem.MolFromSmiles('OCC1CC1'),))

In [12]: Chem.MolToSmiles(ps[0][0])
Out[12]: 'BrCC1CC1'

[thorlucas]

It looks like you are expecting the RDKit reaction code to do something other than what it does.
From the RDKit's perspective that means that each reactant of r2 is a substructure of a reactant of r1 and each product of r2 is a substructure of a product of r1.

It seems to me that it would make sense a blank should match anything. From Daylight's page on SMARTS:

Query:	Target:	Matches:
C>>	CC>>CN	2
>C>	CC>>CN	0
>>C	CC>>CN	1
C	CC>>CN	3

Or is the function in RDKit not meant to mimick Daylight's SMARTS behavior?

That is correct. The substructure matching code currently ignores atom mapping information. This is something we should look into changing, particularly when doing substructure matches on reactions.

Is this a bug or intentional? Looking at the documentation it seems like perhaps internally it loops through all reactants and products and compares them molecule by molecule — and maybe the molecule-molecule query doesn't have atom mapping (becuase that wouldn't make sense), am I right? It seems like you'd have to compare reactions holistically in order for it to work.

[thorlucas]

I'm not sure what you're trying to do.

I'm currently writing a vim plugin that essentially greps all files in a wiki directory, looking for SMARTS/SMILES strings. Then it matches them against a query to see if the reaction is relevant. This is useful because I can ask the questions: "what reactions have I written about that transform X into Y?" or "I need to make Y. What reactions can I use?" or "I have X. What can I do with it?"

It seems like RDKit's focus on SMARTS reactions is to use them for product enumeration. My goal is to use it instead for querying.

So in the case of that specific reaction "[C;D2:1][OH]>>[C;D2:1][Br,Cl,I]" I could have a page on PX3 for example. Then I can ask the question "What takes an alcohol to an alkyl bromide?". More specifically, I could ask "What takes n-butanol to 1-bromobutane? CCCC[OH]>>CCCCBr" and the plugin spits out that page.

But on the flip side, if I were to ask "CCCC[OH]>>CCC(C)CBr", we should not get a match, because the R-groups aren't the same. The reaction string should tell us that the only atom that should change is the oxygen to bromine -- and if the query changes something else, then no match.

I do suppose one roundabout way would be to use what you sent me. For a particular reaction query, I could take the reactants and actually run the reaction, then compare the products. But it seems to me like the functionality I'm describing should already exist, no? Unless I'm mistaking what SMARTS was designed for.

[greglandrum]

It looks like you are expecting the RDKit reaction code to do something other than what it does.
From the RDKit's perspective that means that each reactant of r2 is a substructure of a reactant of r1 and each product of r2 is a substructure of a product of r1.

It seems to me that it would make sense a blank should match anything. From Daylight's page on SMARTS:

Query:
Target:
Matches:

C>>
CC>>CN
2

C>
CC>>CN
0

C
CC>>CN
1

C
CC>>CN
3

Or is the function in RDKit not meant to mimick Daylight's SMARTS behavior?

With the exception of the last example, which doesn't use a reaction as a query, the RDKit behaves EXACTLY like the daylight examples:

In [13]: query = rdChemReactions.ReactionFromSmarts('C>>')

In [14]: rxn = rdChemReactions.ReactionFromSmarts('CC>>CN')

In [15]: rdChemReactions.HasReactionSubstructMatch(rxn,query)
Out[15]: True

In [16]: query2 = rdChemReactions.ReactionFromSmarts('>>C')

In [17]: rdChemReactions.HasReactionSubstructMatch(rxn,query2)
Out[17]: True

you just need to call HasReactionSubstructMatch() with the arguments in the correct order.

In order to avoid further confusion, I'm going to wait until we agree about this basic point before doing anything else on this thread.