Prepare Fragmenstein for bioisosteres

[matteoferla]

To Fix:
* in a two mol constrained mapping: one map seems to be disreguarded.

• Monster and Victor accept a custom_map, but Laboratory does not.

EDIT: there was a mistake in the test

Current issues against 2nd point

1. Monster cannot deal with discontinuous placements, e.g. placing a molecule with a longer linker than in the reference needs the reference to be split into two non-overlapping fragments.
2. Stereoisomers are split into separate entries for inference by `Laboratory.

Aim of 2nd point

@stephwills wants to place a bioisostere according to a reference compound, which would donate pharmacophoric atoms, each non-adjecent atom needs to split into atoms due to 1.

When there are two equal atoms/pharmacophores, eg. oxygen elements in a carboxylate, a bioisostere can be mapped in at least two orientations due to 2.

Major issue lurking for 2nd

In general, the pharmacophores may need better treatment. For now, proving custom_map will do.

Code

This all works

from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Geometry import Point3D
from fragmenstein import Monster
from typing import List

def extract_atoms(mol, zahl: int, prefix: str='pharmacophore') -> List[Chem.Mol]:
    n_query = AllChem.AtomNumEqualsQueryAtom(zahl)
    atom_iter = enumerate(mol.GetAtomsMatchingQuery(n_query))
    return [atom_to_mol(atom, f'{prefix}#{i}') for i, atom in atom_iter]

def atom_to_mol(atom: Chem.Atom, name: str):
    i: int = atom.GetIdx()
    element : str = atom.GetSymbol()
    mol: Chem.Mol = Chem.MolFromSmiles(element)
    point: Point3D = atom.GetOwningMol().GetConformer().GetAtomPosition(i)
    conf = Chem.Conformer(1)
    conf.SetAtomPosition(0,  point)
    mol.AddConformer(conf)
    mol.SetProp('_Name', name)
    return mol
    
# easily paddable
 diaminopentane = AllChem.AddHs(Chem.MolFromSmiles('NCCCCCN'))
AllChem.EmbedMolecule(diaminopentane)
nitros = extract_atoms(diaminopentane, 7)
diaminoethane = Chem.MolFromSmiles('NCCN')
monstah = Monster(nitros)
monstah.journal.handlers = []
monstah.place(diaminoethane, custom_map={'pharmacophore#0': {0: 0}, 'pharmacophore#1': {0: 3}})
assert monstah.convert_origins_to_custom_map()['pharmacophore#0'][0] == 0
assert monstah.convert_origins_to_custom_map()['pharmacophore#1'][0] == 3

glycine = AllChem.AddHs(Chem.MolFromSmiles('OC(=O)CN'))
AllChem.EmbedMolecule(glycine)
oxys = extract_atoms(glycine, 8)
oxaldehyde = Chem.MolFromSmiles('O=CC=O')
monstah = Monster(oxys)
monstah.journal.handlers = []
monstah.place(oxaldehyde, custom_map={'pharmacophore#0': {0: 0}, 'pharmacophore#1': {0: 3}})
print(monstah.convert_origins_to_custom_map())
assert monstah.convert_origins_to_custom_map()['pharmacophore#0'][0] == 0
assert monstah.convert_origins_to_custom_map()['pharmacophore#1'][0] == 3

# real life glycine
glyblock = '''HETATM    1  N   GLY C 306      -1.229   6.136 -32.943  1.00 69.73      S    N  
HETATM    2  CA  GLY C 306      -1.174   7.098 -31.820  1.00 60.45      S    C  
HETATM    3  C   GLY C 306      -2.101   8.290 -32.017  1.00 63.90      S    C  
HETATM    4  O   GLY C 306      -1.680   9.236 -32.745  1.00 59.85      S    O  
HETATM    5  OXT GLY C 306      -3.219   8.254 -31.386  1.00 44.53      S    O  
HETATM    6  H01 GLY C 306      -1.911   5.419 -32.740  1.00 69.73      S    H  
HETATM    7  H02 GLY C 306      -0.155   7.481 -31.766  1.00 60.45      S    H  
HETATM    8  H03 GLY C 306      -1.460   6.581 -30.904  1.00 60.45      S    H  
HETATM    9  H04 GLY C 306      -0.321   5.712 -33.068  1.00 69.73      S    H  
CONECT    1    2    6    9
CONECT    2    1    3    7    8
CONECT    3    2    4    5
CONECT    4    3
CONECT    5    3
CONECT    6    1
CONECT    7    2
CONECT    8    2
CONECT    9    1
END'''
glycine = Chem.MolFromPDBBlock(glyblock)
glycine.SetProp('_Name', 'glycine')
oxys = extract_atoms(glycine, 8)
oxaldehyde = Chem.MolFromSmiles('O=CC=O')
monstah = Monster(oxys)
monstah.journal.handlers = []
monstah.place(acylsulfonimide, custom_map={'pharmacophore#0': {0: 2}, 'pharmacophore#1': {0: 5}})
print(monstah.convert_origins_to_custom_map())
assert monstah.convert_origins_to_custom_map()['pharmacophore#0'][0] == 2
assert monstah.convert_origins_to_custom_map()['pharmacophore#1'][0] == 5

cyclopentanedione = Chem.MolFromSmiles('C1CC(C)=C(-[O-])C1=O')
cyclopentanedione.SetProp('_Name', 'cyclopentanedione')

oxadiazolone = Chem.MolFromSmiles('CC1NC(=O)ON=1')
oxadiazolone.SetProp('_Name', 'oxadiazolone')

acylsulfonimide = Chem.MolFromSmiles('CC(=O)NS(=O)(=O)C')
acylsulfonimide.SetProp('_Name', 'acylsulfonimide')
monstah = Monster(oxys)
monstah.journal.handlers = []
monstah.place(acylsulfonimide, custom_map={'pharmacophore#0': {0: 2}, 'pharmacophore#1': {0: 5}})
print(monstah.convert_origins_to_custom_map())
assert monstah.convert_origins_to_custom_map()['pharmacophore#0'][0] == 2
assert monstah.convert_origins_to_custom_map()['pharmacophore#1'][0] == 5