Mapping and merging issues

lohedges commented 5 years ago

This issue is to document mapping and merging issues for the Cathepsin ligands from the D3R2018 data set.

The table below outlines the suggested ligand pairs to be used in the "soft-go" BioSimSpace trial. The table indicates the name of the ligand pairs from the test set (numbers in brackets are the number of atoms in the ligand), the number of atom pair matches found by BioSimSpace using its default settings, and whether or not it was possible to merge the two molecules based on the mapping.

Ligand 1	Ligand2	Matches	Merges
CatS_178 (94)	CatS_29 (96)	45	No
CatS_29 (96)	CatS_178 (94)	45	No
CatS_178 (94)	CatS_191 (100)	58	No
CatS_191 (100)	CatS_178 (94)	58	No
CatS_191 (100)	CatS_30 (92)	63	Yes
CatS_30 (92)	CatS_191 (100)	63	Yes
CatS_30 (92)	CatS_29 (96)	89	Yes
CatS_29 (96)	CatS_30 (92)	89	Yes
CatS_155 (93)	CatS_30 (92)	92	Yes
CatS_30 (92)	CatS_155 (93)	92	Yes
CatS_155 (93)	CatS_153 (90)	86	Yes
CatS_153 (90)	CatS_155 (93)	86	Yes
CatS_153 (90)	CatS_157 (92)	90	Yes
CatS_157 (92)	CatS_153 (90)	90	Yes
CatS_157 (92)	CatS_155 (93)	86	Yes
CatS_155 (93)	CatS_157 (92)	86	Yes

As seen, there are several ligand pairs that cannot be merged successfully. These are the ones with the lowest number of matched atoms. In all cases, the maximum common substructure (MCS) mapping is symmetric with respect to the order of the ligands.

To reproduce any of the above data, navigate to CatS/BSS/ligands_aligned/parameterised and run the following (modifying the ligand names accordingly for the pair of interest):

import BioSimSpace as BSS

# Load the two ligands.
m0 = BSS.IO.readMolecules(BSS.IO.glob("para*/CatS_155.*")).getMolecules()[0]
m1 = BSS.IO.readMolecules(BSS.IO.glob("para*/CatS_157.*")).getMolecules()[0]

# Generate the forward and reverse mappings.
mapping0 = BSS.Align.matchAtoms(m0, m1)
mapping1 = BSS.Align.matchAtoms(m1, m0)

# Assert the mappings are the same length.
assert len(mapping0) == len(mapping1)

# Align the molecules based on the mapping.
m0_aligned = BSS.Align.rmsdAlign(m0, m1, mapping0)
m1_aligned = BSS.Align.rmsdAlign(m1, m0, mapping1)

# Now try to merge the molecules based on the mappings.
try:
    merged0 = BSS.Align.merge(m0_aligned, m1, mapping0)
except:
    raise

try:
    merged1 = BSS.Align.merge(m1_aligned, m0, mapping1)

    # Write the merged molecule to PDB format.
    BSS.IO.saveMolecules("merged1", merged1, "pdb")
except:
    raise

# Write the lambda = 0 state of the merged molecules to PDB format.

map0 = {}
for prop in merged0._sire_molecule.propertyKeys():
    if prop[-1] == "0":
        map0[prop[:-1]] = prop
 BSS.IO.saveMolecules("merged0", merged0, "pdb", property_map=map0)

map1 = {}
for prop in merged1._sire_molecule.propertyKeys():
    if prop[-1] == "0":
        map1[prop[:-1]] = prop
 BSS.IO.saveMolecules("merged1", merged1, "pdb", property_map=map1)

lohedges commented 5 years ago

If anyone has any intuition as to why the ones that won't merge are failing, that would be great. In my previous experience I come across merging issues when there are mismatches in ring sizes, or partial ring overlaps following the alignment. Since our alignment scheme only performs simple translations and rotations it's possible that the resulting structure is partially misformed, or has some clashes.

jmichel80 commented 5 years ago

Hi Lester,

Could you update the above script such that merge molecules are written down as pdb or mol2 files for ease of visualisation ? I need to see structures to use my 'intuition' :-) This will make it easier to check that the cases where merged did not raise an exception are actually valid merged molecules. In particular this pair is suspicious

CatS_191 (100) CatS_30 (92) 63 Yes CatS_30 (92) CatS_191 (100) 63 Yes

BW

Julien

Dr. Julien Michel, Senior Lecturer Room 263, School of Chemistry University of Edinburgh David Brewster road Edinburgh, EH9 3FJ United Kingdom phone: +44 (0)131 650 4797 http://www.julienmichel.net/

On Thu, Oct 25, 2018 at 3:00 PM Lester Hedges notifications@github.com<mailto:notifications@github.com> wrote:

If anyone has any intuition as to why the ones that won't merge are failing, that would be great. In my previous experience I come across merging issues when there are mismatches in ring sizes, or partial ring overlaps following the alignment. Since our alignment scheme only performs simple translations and rotations it's possible that the resulting structure is partially misformed, or has some clashes.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://github.com/michellab/D3R2018/issues/1#issuecomment-433062093, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ALLd5BVbFtG-ju2JVluaPuKh3hvLoE9Dks5uocO-gaJpZM4X6VEd.

The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.

lohedges commented 5 years ago

The script has been updated. I just make use of a property map to write the lambda = 0 state of each merged molecule. The same thing can be done to output the lambda = 1 state too. (Just replace if prop[-1] == "0" with if prop[-1] == "1".)

lohedges commented 5 years ago

I also think Toni has done a good job of pre-aligning the ligands. You should just be able to navigate to CatS/BSS/ligands_aligned and load any two of PDB files corresponding the pairs above in order to see what the alignment looks like.

Alternatively, you can save the aligned ligands generated by BioSImSpace and compare those, e.g.

BSS.IO.saveMolecules("aligned0", m0_aligned, "pdb")
BSS.IO.saveMolecules("aligned1", m1_aligned, "pdb")

michellab / D3R2018

Mapping and merging issues #1

Dr. Julien Michel, Senior Lecturer Room 263, School of Chemistry University of Edinburgh David Brewster road Edinburgh, EH9 3FJ United Kingdom phone: +44 (0)131 650 4797 http://www.julienmichel.net/