lohedges
commented
1 year ago Thanks, I'll take a closer look when I get a chance. We had quite a few tests for weird ring-size and ring-breaking edge cases, but unfortunately most of these now need to be marked as xfail due to changes in the RDKit API generating different mappings. If I work out what's going on here, then I'll add these systems as examples. It's quite hard to figure out all of the possible ways in which this can/can't occur and essentially involves manually looking at the connectivity of the molecule (as based on the mapping) and working things out.
Out of interest, do you have the original SDF files used to load the ligands, i.e. prior to parameterisation? I wonder whether these would give the same mappings. (The SDF specific attributes would be retained as properties of the molecules.) We now use Sire's built in to_rdkit() functionality when performing the MCS, but this might be more robust if going via a molecule originally loaded from SDF, e.g. if stererochemistry is important, or can't be inferred any other way.
Just to note that when converting to RDKit with Sire (say for ligands 27 and 42) I get valence warnings, and if I attempt to compute a mapping using the LomaAtomMapper from openfe then I get the same warning and no mapping.
In [1]: import sire as sr
In [2]: from gufe import SmallMoleculeComponent
In [3]: from openfe.setup.atom_mapping import LomapAtomMapper
In [4]: m0 = sr.load("ligands/*27*")
In [5]: m1 = sr.load("ligands/*42*")
In [6]: rdm0 = sr.convert.to_rdkit(m0)
[15:10:12] Explicit valence for atom # 1 C, 5, is greater than permitted
In [7]: rdm1 = sr.convert.to_rdkit(m1)
[15:10:14] Explicit valence for atom # 8 C, 5, is greater than permitted
In [8]: openfe_m0 = SmallMoleculeComponent(rdm0)
In [9]: openfe_m1 = SmallMoleculeComponent(rdm1)
In [10]: mapper = LomapAtomMapper()
In [11]: mapping_gen = mapper.suggest_mappings(openfe_m0, openfe_m1)
In [12]: mapping = next(mapping_gen)
[15:10:53] Explicit valence for atom # 1 C, 5, is greater than permitted
╭─────────────────────────────── Traceback (most recent call last) ────────────────────────────────╮
│ in <module>:1 │
╰──────────────────────────────────────────────────────────────────────────────────────────────────╯
StopIteration
annamherz
Hello! I'm currently matching and merging molecules for MCL1 ligands. The basic code used to replicate these is:
With the files named as the ligands attached. ligands.zip
There are some different situations that occur in relation to ring systems, outlined below.
Case 1
In the most basic case, as in the case of some ligands (lig_27~lig_42), this does not proceed correctly initially. Many of the atoms are included in the perturbable region and some of the shared atoms are not included in the MCS.
When the mapping dictionary
Notably, this is able to proceed with
{21:6}is added asprematch, this is then able to proceed okay. Prematch for this ligand series is chosen so the N in the core region match.ring_size=Falseandring_break=Falseand the ring is mapped correctly to the fused ring system, which is not the case in Case 2.Case 2
However in a another case, lig_27~lig_43, when the mapping dictionary
I think ideally, with the introduced mapping, this should still be able to proceed with ring_size and ring_break set as False, and the whole ring region in this case is taken as the perturbed region in each case. So for lig_27 the phenyl ring and for lig_43 the entire fused ring system. Then if these were both set to True, in that case the phenyl ring should ideally be mapped onto an entire ring in the fused phenyl as opposed to split between the two fused rings.
The difference between lig_42 and lig_43 in this case is minor, with lig_42 being a 1-napthyl group and lig_43 being a 2-napthyl group, but this causes a large difference in the mapping.
{ 21:26 }which is also the N is added, this is not able to proceed without settingring_size=True,ring_break=True. In this case then the resulting perturbed region is like below:Case 3
Here, lig_43~lig_45, there is a change of hybridisation state with the rings. Again, without prematch, the core region of the molecule is ignored:
And when prematch
{26:6}is introduced, withring_size=False,ring_break=Falseit then maps the differently hybridised C to each other. The differently hybridised atoms should not map to each other as this is not great for the free energy perturbation.I was wondering if there would be some way to programmatically address the mapping and merging ring perturbations so the intended perturbable molecule can be obtained? I guess the atoms could also be removed manually from the generated mapping that is passed into the merge function, however this requires quite some time per perturbation to do manually as each atom needs to be inspected then.
Thank you!