Questions about mapping

[SZhang66]

Hi,

Recently, I've been trying to set up the RBFE calculations for the benchmark systems cmet and eg5 using perses (NES). I'd like to see if we could use it to study the transformation between charged compounds. But I encountered two issues, and I think maybe you know the reason since you're more experienced. Appreciate for your time and any help!

1.) For the edge 15 -> 22 of cmet, I got below A-mapping when I run perses.app.relative_setup.RelativeFEPSetup(...) to setup the system. But if I run perses.rjmc.atom_mapping.AtomMapper(), I could get B-mappinng which looks like more reasonable. As I understand it, RelativeFEPSetup is also using rjmc.atom_mapping to create a mapping. So, I have no idea why they're so different.

• FYI, 15 and 22 are the indices of ligands in the attached inputs/cmet_ligands.sdf. And the script (inputs/check_mappings.py)I run for quick testing is also attached.

• A-mapping

• B-mapping

2.) The second issue is about the edge 1->7 of eg5: I got below error msg when I run perses.app.relative_setup.RelativeFEPSetup(...) to setup the system. As long as the transformation includes ligand 1 or ligand 7, I can see the same error. My guess is that the pyramidal nitrogen connecting different R1,R2,R3 in the 1st and 7th ligands is causing this issue. But I'm not 100% sure.

• FYI, the script (inputs/check_mappings.py) is also used for this test.

  openff.toolkit.utils.exceptions.UndefinedStereochemistryError: Unable to make OFFMol from OEMol: OEMol has unspecified stereochemistry. oemol.GetTitle(): MOL
  Problematic atoms are:
  Atom atomic num: 7, name: N2, idx: 3, aromatic: False, chiral: True with bonds:
  bond order: 1, chiral: False to atom atomic num: 6, name: C2, idx: 2, aromatic: False, chiral: False
  bond order: 1, chiral: False to atom atomic num: 6, name: C3, idx: 4, aromatic: False, chiral: False
  bond order: 1, chiral: False to atom atomic num: 6, name: C23, idx: 29, aromatic: True, chiral: False

Attached files: inputs.zip outputs.zip Pkgs used: perses 0.10.1
openff-toolkit 0.10.0

[ijpulidos]

Thanks for the feedback. I can reproduce the issue, I'm going to dig deeper to see what may be causing it, because it seems like an important bug. It seems that there's something happening along the way to the ligands before the topology proposal is being made.

[ijpulidos]

@SZhang66 The problem seems to happen with using OEMol objects from openeye-toolkits. I thought updating the version to 2022.2.2 would solve it but that wasn't the case. You can test the following code to see that mapping using the openeye molecules gives the suboptimal mapping. I see this issue with versions 2021.2.0 and 2022.2.2 of the openeye-toolkits

from openff.toolkit.topology import Molecule
from perses.rjmc.atom_mapping import AtomMapper

def mapping(ligand_input, old_ligand_index, new_ligand_index):
    mols = Molecule.from_file(ligand_input)
    old_mol = mols[old_ligand_index]
    new_mol = mols[new_ligand_index]

    atom_mapper = AtomMapper()
    atom_mapping = atom_mapper.get_best_mapping(old_mol, new_mol)
    score = atom_mapper.score_mapping(atom_mapping) 
    print ('score=', score)
    atom_mapping.render_image('best_mapping_offmols.png')    # save mapping image

def mapping_oemols(ligand_input, old_ligand_index, new_ligand_index):
    from openeye import oechem
    ifs = oechem.oemolistream(ligand_input)
    mols_list = [oechem.OEMol(mol) for mol in ifs.GetOEMols()]
    old_mol = mols_list[old_ligand_index]
    new_mol = mols_list[new_ligand_index]

    # Do mapping before relative setup with oe mols
    atom_mapper = AtomMapper()
    atom_mapping = atom_mapper.get_best_mapping(old_mol, new_mol)
    score = atom_mapper.score_mapping(atom_mapping)
    print('score=', score)
    atom_mapping.render_image("best_mapping_oemols.png")

#------
ligand_input = "cmet_ligands.sdf"   #eg5_ligands.sdf"
old_ligand_index = 15    #1
new_ligand_index = 22    #7

mapping(ligand_input, old_ligand_index, new_ligand_index)
mapping_oemols(ligand_input, old_ligand_index, new_ligand_index)

[SZhang66]

@ijpulidos Thanks for testing it! And Yes, I got the exact same problems when I run the above code you provided using openeye-toolkits 2021.2.0and 2022.1.1. Guess some parameters can't be recognized or wrong atoms are identified for these compounds? Really curious what caused it?

[mikemhenry]

Maybe related to #1161

[ijpulidos]

Good news, as far as I can see this has been solved by the changes in #1128 :tada: which is already merged in main

[SZhang66]

Tested with the latest update. Now it works fine with these edges. Thanks!