multiple smarts in --cut-smarts #15
Comments
|
Andrew commented on this request on the RDKit-discuss mailing list as: The SMARTS handling code only touches <50 lines of code. It does not seem that hard to have it take multiple --cut-smarts, apply each of the cuts, find the unique union of those cuts, and work with them. Could you add that as a issue in the mmpdb tracker? It is in principle possible to merge two fragment files together and index the result. However, it would be difficult to use the indexed database for analysis purposes, because any input/query structure would use the single SMARTS pattern defined in the database. |
|
At the RDKit UGM Hackathon 2019, this question came up again. Participants wanted to use the RECAP rules for cutting. Creating a single SMARTS to match all 11 rules might theoretically be possible, but would results in an extremely complicated string which would then be hard to debug and modify. Extending RDKit such that a list of SMARTS is appears as the preferred long term solution. |
|
There are 12 rules, not 11: How to people want to specify the cut with these? Is the cut match defined with the product side of the reaction, and the reactant side ignored? Some of those SMARTS use more than two atoms. The first makes a cut between :1 and :2 while the second makes a cut between :2 and :3. That means that if the reaction side is ignored (eg, if the cut is always made between :1 and :2) then there will be problems. It could do a more in-depth analysis of the transform to detect if there is a labeled pair on the product side which is not a labeled pair in the reactant side, and use that for the cut. But that's overkill if people really just want |
|
I'm thinking to support it as Looking at the RECAP rules, there are several places where I see problems.
Given NC(=O)N this removes the C(=O) to give N.N. Should the SMARTS be
The existing code only allows cuts on single bonds. The RECAP patterns use Note that There are far more matches with It looks like in those few cases where the non-ring bond type is not specified, it's okay for me to say it's a single bond, without changing the intent of matching an amide. It also seems like that RECAP definition in RDKit is wrong, in that it is not supposed to match a double bond there.
The pattern mmpdb cannot handle this case. Should I drop it? |
|
Going back to acquaregia's request, can you give an example of of the SMIRKS you are interested in? I can see two steps that might be affected: 1) limit the fragment to just a few SMARTS patterns, and 2) limit the indexing to just a few SMIRKS patterns. I would like to see some of the SMIRKS to get a better feel for how to handle this. For example, if all of the SMIRKS were transforms of R-groups to R-groups, where the R-groups could be defined as SMILES fragments with a single attachment point denoted Otherwise, if multiple distinct SMARTS are needed, then the mmpdb file formats need to change someone in order to store them. There could be multiple entries, one per definition, or they could be space/tab separated. |
|
@adalke, I know you've stepped back from this, so this is more a record of my thoughts than a request... I came across this issue today, when trying to work out why the database was 'missing' a matched pair compared to what I expected, and it led me to consider how the available fragmentation patterns differ from how Matsy worked. The Matsy fragmentation scheme is described in https://pubs.acs.org/doi/10.1021/jm500022q:
This came from Roger. It's deceptively simple, but captures the the synthetically aware sense of the RECAP rules. Subsequently, or at the same time, Antonio and Bajorath did work in this area explicitly using RECAP (here's one reference, https://pubs.rsc.org/en/content/articlelanding/2014/md/c3md00259d). We didn't use SMARTS, but I think the following is equivalent in mmpdb terms: I've been scratching my head wondering how to use the existing codebase with these patterns. Can I combine these into a single gnarly recursive SMARTS, where each end of the bond recursively contains the full pattern? It's not going to be pretty and it's not going to be fast - that's for sure. Since I'm only interested in indexing (i.e. identifying matched pairs), it sounds like I can build two separate dbs and merge them. I'll try that first. |
|
On Sep 5, 2024, at 11:16, baoilleach ***@***.***> wrote:
I've been scratching my head wondering how to use the existing codebase with these patterns. Can I combine these into a single gnarly recursive SMARTS, where each end of the bond recursively contains the full pattern? It's not going to be pretty and it's not going to be fast - that's for sure. Since I'm only interested in indexing (i.e. identifying matched pairs), it sounds like I can build two separate dbs and merge them. I'll try that first.
They cannot be combined now.
One option is to create a superset SMARTS grammar, eg, use "%%" to separate SMARTS patterns ("%" cannot be used to start or end a SMARTS, and "%%" is not valid in SMARTS; alternative options exist, like the space character)
- Change fragment_types.py in parse_cut_smarts() so it returns a list of compiled SMARTS molecules, rather than a single.
smarts_terms = smarts.split("%%")
if not smarts_terms:
raise ValueError("cut smarts must not be empty")
patterns = []
for smarts_term in smarts_terms:
pattern = Chem.MolFromSmarts(smarts_term)
if pattern is None:
raise ValueError("unable to parse cut SMARTS")
if pattern.GetNumAtoms() != 2:
raise ValueError("cut SMARTS must match exactly two atoms")
if pattern.GetNumBonds() != 1:
raise ValueError("cut SMARTS must connect both atoms")
patterns.append(pattern)
return patterns
a fancier one to give a more narrowed down error reporting might be:
smarts_terms = smarts.split("%%")
if not smarts_terms:
raise ValueError("cut smarts must not be empty")
def get_value_error(msg, smarts_term):
# if there are multiple terms then narrow it down
if len(smarts_terms) == 1:
extra = ""
else:
extra = " term {smarts_term!r}"
return ValueError(msg.format(extra=extra))
patterns = []
for smarts_term in smarts_terms:
pattern = Chem.MolFromSmarts(smarts_term)
if pattern is None:
raise get_value_error("unable to parse cut SMARTS{extra}", smarts_term)
if pattern.GetNumAtoms() != 2:
raise get_value_error("cut SMARTS{extra} must match exactly two atoms", smarts_term)
if pattern.GetNumBonds() != 1:
raise get_value_error("cut SMARTS{extra} must connect both atoms", smarts_term)
patterns.append(pattern)
return patterns
- Change fragment_types.py "cut_pattern" to "cut_patterns" to reflect the change.
- Change fragment_types.py "get_cut_atom_pairs" to handle the new logic
seen = set()
for pat in self.cut_patterns:
for (atom1_idx, atom2_idx) in mol.GetSubstructMatches(self.cut_pattern):
# put into canonical order so cuts are consistent across all patterns
if atom1_idx < atom2_idx:
seen.add((atom1_idx, atom2_idx))
else:
seen.add((atom2_idx, atom1_idx))
return list(seen)
This would, I think, allow --cut-smarts on the command-line, and stored in the database, and be backwards compatible.
Andrew
***@***.***
|
|
Works like a dream:: I've checked it in over at https://github.com/baoilleach/mmpdb. |
|
On Sep 5, 2024, at 17:35, baoilleach ***@***.***> wrote:
Works like a dream::
Good to hear! I like how it works to extend existing names.
In retrospect, I think "||" is a better separator than "%%" as "|" is not used in SMARTS at all, as "|" has a meaning as "or" , and "||" has an even stronger meaning as "or".
Andrew
***@***.***
|
|
Indeed. Done. |
I have large database 500K compounds and I am interested in finding only few transforms.
Ideally I would like to give transform in the form of smirks.
I understand that it might be easier to ask for a different fragmentation pattern and perform indexing on it.
I can translate the smirks into smarts specifying specific bonds.
For the tool to be useful I would like to be able to provide more than one SMARTS to the --cut-smarts option.
It would be excellent if an option like --cache would allow using a fragmentation file and enhance it by specifying other cut patterns.
Thanks.
marco
The text was updated successfully, but these errors were encountered: