[WIP] Biopolymer rdkit smarts

[richardjgowers]

FAO @j-wags

This is how I'd propose doing the standard residue recognition with a cheminformatics toolkit. It's not 100% working yet, I think I might have to tweak the smarts patterns, but this might also improve them (e.g. could do something like [#7;D3;H1] to get the non-cap N ends, and this might resolve the race condition on patterns that there's currently a TODO about). I also think my attempt at string mashing to get the smarts to generalise is weak and should be replaced with some TOOLKIT_REGISTRY.fuzzy_match(SMARTS) to generate queries that are bond/charge/arom insensitive, will have to check OEChem support for that first.

Some questions:

• Before I dive into hacking the smarts patterns, is this approach viable?
• Is the general future outlook for polymer-likes that people can provide a list of monomer-SMARTS to load molecules and apply sufficient information to them?
• Out of curiosity is anyone working on general polymers yet?

[j-wags]

Thanks for putting this together @richardjgowers!

TL;DR

Before I dive into hacking the smarts patterns, is this approach viable?

Yes, I think this is viable. The big pluses are:

• This may be significantly faster (though I can't get this to run locally on T4, so I'd like to verify this!)
• It seems cleaner to use a cheminformatics toolkit for molecule stuff
• It removes the jank around negative atomic nums for hydrogen
• It may allow us to directly query on number of hydrogens/degree

The functional liabilities with the current code, which I think can be totally fixed, are:

• The current code seems to stray away from handling inter-monomer bonds. We do eventually intend to support general polymers and one day the inter-monomer bond may not be single. I don't think this is a fundamental limitation of the approach, though, and it shouldn't be too hard to match these (if it's not happening already)
• The removal of checks for unassigned edges makes me nervous - Even if we let peptide bonds slide, we may be missing a sidechain and not knowing it

Details

I might have to tweak the smarts patterns, but this might also improve them (e.g. could do something like [#7;D3;H1] to get the non-cap N ends, and this might resolve the race condition on patterns that there's currently a TODO about).

Oh, I really like that we could get rid of that race condition, and start explicitly stating the number of hydrogens and atom degrees.

The code that generates the substructure dict should be generalizable enough to write SMARTS with that info, though the exact implementation may take some thought. I'll need to stew on this a little bit but could probably get that working if we choose to pursue this.

I also think my attempt at string mashing to get the smarts to generalise is weak and should be replaced with some TOOLKIT_REGISTRY.fuzzy_match(SMARTS) to generate queries that are bond/charge/arom insensitive, will have to check OEChem support for that first.

Yup, pretty much agree with all of this.

Thinking a bit more about the method signature above - It should probably be TOOLKIT_REGISTRY.fuzzy_match(SMARTS, molecule), but where the molecule isn't yet chemically valid. This is tricky because our Molecule class isn't made to hold chemically invalid molecules, so this could get brittle (hence why I was using a networkx graph). Maybe the method could take an entire substructure dict and an OpenMM topology, and then let the toolkit wrapper do the complete assignment and return an OFFMol?

Is the general future outlook for polymer-likes that people can provide a list of monomer-SMARTS to load molecules and apply sufficient information to them?
Out of curiosity is anyone working on general polymers yet?

Yes. @connordavel is an undergrad in the Shirts group currently looking at handling more complex polymers. Right now the "substructure dictionary" files are "private" (we don't guarantee their stability or format) and Connor is hacking around directly on the contents. We will make them public eventually, but not in the 0.11.0 release (case in point, this PR might rejigger the substructure dict to query on n_hydrogens).

[richardjgowers]

@j-wags I've updated this branch to be my best attempt so far, my working gist here: https://gist.github.com/richardjgowers/c40186e1fe18ebafb22958e18827dd2f

It's getting most but not all AAs, using the patterns in the existing json (and eventually we'll want a synthetic test that iterates through the json and creates one of each monomer to check we correctly assign to it). If I change the input patterns to not include explicit hydrogens it seems to work, but then that's annoying if you want to assign atom names to the hydrogens as they won't show up in the resulting match.

[j-wags]

Thanks for documenting your progress, @richardjgowers - I'll take over this PR!

[j-wags]

I'm going to repurpose this PR to port these changes into a branch rather than a fork so CI will be happy. I'll open a new PR momentarily!