drrdom / crem Goto Github PK

View Code? Open in Web Editor NEW

193.0 7.0 37.0 427 KB

CReM: chemically reasonable mutations framework

License: BSD 3-Clause "New" or "Revised" License

Python 3.54% Shell 0.03% Jupyter Notebook 96.43%

chemoinformatics chemical-space

crem's People

Contributors

Stargazers

Watchers

crem's Issues

Multiple replacements in grow_mol

Hi, thanks for creating this useful repo which allows us to generate new molecules from a framework!

I have been playing around with the grow_mol function, and I want to specify 3 positions in my starting molecule to grow the molecule. I can do so with the parameter replace_ids and specify the 3 positions in a list. However, I realised that all of the generated molecules only had 1 replacement per molecule. What I am looking for is a way to possibly generate new molecules with 2 or even all 3 positions replaced at once.

I have read the documentation for crem, and as far as I'm aware, there is no way to do this currently. Do you have any advice regarding this issue?

Thank you!
Alan

Parameters relevant for speed in mutate_mol

Hi! Thanks for maintaining this repo! I was wondering which parameters were relevant for speeding up mutate_mol? I've found better speeds by increasing n_cores but if I wanted to speed it up even more would decreasing radius help too?

Database Visualization and Variation Inquiry

Thank you for creating this valuable tool. I'm wondering how I can visualize the database I've created following the documentation. Also, could you please summarize the variations I can expect in the outcome database if I've created it with different radii (for example: 1, 2, 3, 4, 5)? Additionally, what if I only want to create fragmented databases with different heteroatoms and not couple them?

No such table: radius3

Hi! I followed the instruction to install the crem package. However, when I ran the demo code, I found the below issues:

---------------------------------------------------------------------------
OperationalError                          Traceback (most recent call last)
/tmp/ipykernel_1568311/2420725268.py in <module>
      1 m = Chem.MolFromSmiles('c1cc(OC)ccc1C')  # methoxytoluene
----> 2 mols = list(mutate_mol(m, db_name='replacements.db', max_size=1,radius=3))

/home/shuangjia/anaconda3/lib/python3.7/site-packages/crem/crem.py in mutate_mol(mol, db_name, radius, min_size, max_size, min_rel_size, max_rel_size, min_inc, max_inc, max_replacements, replace_cycles, replace_ids, protected_ids, symmetry_fixes, min_freq, return_rxn, return_rxn_freq, return_mol, ncores, **kwargs)
    493                                                                 protected_ids_1=protected_ids, protected_ids_2=None,
    494                                                                 min_freq=min_freq, symmetry_fixes=symmetry_fixes,
--> 495                                                                 **kwargs):
    496             for smi, m, rxn in __frag_replace(mol, None, frag_sma, core_sma, radius, ids, None):
    497                 if max_replacements is None or len(products) < (max_replacements + 1):  # +1 because we added source mol to output smiles

/home/shuangjia/anaconda3/lib/python3.7/site-packages/crem/crem.py in __gen_replacements(mol1, mol2, db_name, radius, dist, min_size, max_size, min_rel_size, max_rel_size, min_inc, max_inc, max_replacements, replace_cycles, protected_ids_1, protected_ids_2, min_freq, symmetry_fixes, **kwargs)
    342                 max_atoms = num_heavy_atoms + max_inc
    343 
--> 344                 row_ids = __get_replacements_rowids(cur, env, dist, min_atoms, max_atoms, radius, min_freq, **kwargs)
    345 
    346                 if max_replacements is None:

/home/shuangjia/anaconda3/lib/python3.7/site-packages/crem/crem.py in __get_replacements_rowids(db_cur, env, dist, min_atoms, max_atoms, radius, min_freq, **kwargs)
    284         else:
    285             sql += f" AND {k} = {v}"
--> 286     db_cur.execute(sql)
    287     return set(i[0] for i in db_cur.fetchall())
    288 

OperationalError: no such table: radius3

I'm not sure if it's a bug after the package update, whether I use pypi to install or manually install this problem exists, look forward to your answer, thank you!

Reproducible output

Implement random seed to make results reproducible

Fragmentation

Was a bit confused by the fragment_mol method in fragmentation.py, specifically these lines:

frags = rdMMPA.FragmentMol(mol, pattern="[!#1]!@!=!#[!#1]", maxCuts=4, resultsAsMols=False, maxCutBonds=30)
frags += rdMMPA.FragmentMol(mol, pattern="[!#1]!@!=!#[!#1]", maxCuts=3, resultsAsMols=False, maxCutBonds=30)

Shouldn't maxCuts=3 be included as a subset of all fragmentations for maxCuts=4?
What exactly is the difference between maxCutBonds and maxCuts? I tried looking for documentation on this and was unfortunately unable to find it.
Would it make more sense for the maxCuts parameter to be user-defined rather than preset?

Thanks for your help!

grow_mol error when used on fragments with dummy atoms

Hello,

Thank you for maintain CReM! I ran into an error while using the program, and would like to see if you have any insights into this problem.

I'm trying to grow a fragment generated using RGroupDecomposition in rdkit, like so:

from rdkit import Chem
from rdkit.Chem import rdRGroupDecomposition as rgd
from crem.crem import grow_mol

mol = Chem.MolFromSmiles("Cc4cccc(Cc3cnc(NC(=O)c1cnn2ccc(C)nc12)s3)c4")
core = Chem.MolFromSmiles("Cc1cnc(NC(=O))s1")
params = rgd.RGroupDecompositionParameters()
params.RGroupLabels = rgd.RGroupLabels.AtomIndexLabels
fragment = rgd.RGroupDecompose([core], [mol], options=params)[0][0]['R1']
display(fragment)

The R1:1 is a dummy atom generated by RGroupDecomposition. If I try to run grow_mol at this point, I receive this error:

derivatives = list(grow_mol(r, db_name='/home/javi/data/rgroup/replacements02_sa2.db', max_replacements=9, return_mol=True))
>>>
ValueError                                Traceback (most recent call last)
crem.ipynb Cell 3' in <cell line: 62>()
---> 66 derivatives = list(grow_mol(r, db_name='replacements02_sa2.db', max_replacements=9, return_mol=True))

File ~/workspace/sci-dev/sci-dev-venv/lib/python3.8/site-packages/crem/crem.py:487, in mutate_mol(mol, db_name, radius, min_size, max_size, min_rel_size, max_rel_size, min_inc, max_inc, max_replacements, replace_cycles, replace_ids, protected_ids, symmetry_fixes, min_freq, return_rxn, return_rxn_freq, return_mol, ncores, **kwargs)
    483 protected_ids = sorted(protected_ids)
    485 if ncores == 1:
--> 487     for frag_sma, core_sma, freq, ids in __gen_replacements(mol1=mol, mol2=None, db_name=db_name, radius=radius,
    488                                                             min_size=min_size, max_size=max_size,
    489                                                             min_rel_size=min_rel_size, max_rel_size=max_rel_size,
    490                                                             min_inc=min_inc, max_inc=max_inc,
    491                                                             max_replacements=max_replacements,
    492                                                             replace_cycles=replace_cycles,
    493                                                             protected_ids_1=protected_ids, protected_ids_2=None,
    494                                                             min_freq=min_freq, symmetry_fixes=symmetry_fixes,
    495                                                             **kwargs):
    496         for smi, m, rxn in __frag_replace(mol, None, frag_sma, core_sma, radius, ids, None):
    497             if max_replacements is None or (max_replacements is not None and len(products) < max_replacements):

File ~/workspace/sci-dev/sci-dev-venv/lib/python3.8/site-packages/crem/crem.py:339, in __gen_replacements(mol1, mol2, db_name, radius, dist, min_size, max_size, min_rel_size, max_rel_size, min_inc, max_inc, max_replacements, replace_cycles, protected_ids_1, protected_ids_2, min_freq, symmetry_fixes, **kwargs)
    334 hac_ratio = num_heavy_atoms / mol_hac
    336 if (min_size <= num_heavy_atoms <= max_size and min_rel_size <= hac_ratio <= max_rel_size) \
    337         or (replace_cycles and cycle_pattern.search(core)):
--> 339     frag_sma = combine_core_env_to_rxn_smarts(core, env)
    341     min_atoms = num_heavy_atoms + min_inc
    342     max_atoms = num_heavy_atoms + max_inc

File ~/workspace/sci-dev/sci-dev-venv/lib/python3.8/site-packages/crem/mol_context.py:349, in combine_core_env_to_rxn_smarts(core, env, keep_h)
    346         links[i].append(a.GetNeighbors()[0].GetIdx())
    347         att_to_remove.append(a.GetIdx())
--> 349 for i, j in links.values():
    350     m.AddBond(i, j, Chem.BondType.SINGLE)
    352 for i in sorted(att_to_remove, reverse=True):

ValueError: too many values to unpack (expected 2)

However, when I try to replicate the problem later on in a notebook, I get a different error:

list(grow_mol(fragment, db_name='replacements02_sa2.db', max_replacements=9, return_mol=True))
>>>
RuntimeError                              Traceback (most recent call last)
crem.ipynb Cell 4' in <cell line: 10>()
---> 10 list(grow_mol(fragment, db_name='replacements02_sa2.db', max_replacements=9, protected_ids=[9,10], return_mol=True))

File ~/workspace/sci-dev/sci-dev-venv/lib/python3.8/site-packages/crem/crem.py:487, in mutate_mol(mol, db_name, radius, min_size, max_size, min_rel_size, max_rel_size, min_inc, max_inc, max_replacements, replace_cycles, replace_ids, protected_ids, symmetry_fixes, min_freq, return_rxn, return_rxn_freq, return_mol, ncores, **kwargs)
    483 protected_ids = sorted(protected_ids)
    485 if ncores == 1:
--> 487     for frag_sma, core_sma, freq, ids in __gen_replacements(mol1=mol, mol2=None, db_name=db_name, radius=radius,
    488                                                             min_size=min_size, max_size=max_size,
    489                                                             min_rel_size=min_rel_size, max_rel_size=max_rel_size,
    490                                                             min_inc=min_inc, max_inc=max_inc,
    491                                                             max_replacements=max_replacements,
    492                                                             replace_cycles=replace_cycles,
    493                                                             protected_ids_1=protected_ids, protected_ids_2=None,
    494                                                             min_freq=min_freq, symmetry_fixes=symmetry_fixes,
    495                                                             **kwargs):
    496         for smi, m, rxn in __frag_replace(mol, None, frag_sma, core_sma, radius, ids, None):
    497             if max_replacements is None or (max_replacements is not None and len(products) < max_replacements):

File ~/workspace/sci-dev/sci-dev-venv/lib/python3.8/site-packages/crem/crem.py:314, in __gen_replacements(mol1, mol2, db_name, radius, dist, min_size, max_size, min_rel_size, max_rel_size, min_inc, max_inc, max_replacements, replace_cycles, protected_ids_1, protected_ids_2, min_freq, symmetry_fixes, **kwargs)
    312 else:
    313     mol = mol1
--> 314     f = __fragment_mol(mol, radius, protected_ids=protected_ids_1, symmetry_fixes=symmetry_fixes)
    316 if f:
    317     mol_hac = mol.GetNumHeavyAtoms()

File ~/workspace/sci-dev/sci-dev-venv/lib/python3.8/site-packages/crem/crem.py:95, in __fragment_mol(mol, radius, return_ids, keep_stereo, protected_ids, symmetry_fixes)
     92         atom.SetIntProp("Index", atom.GetIdx())
     94 # heavy atoms
---> 95 frags = rdMMPA.FragmentMol(mol, pattern="[!#1]!@!=!#[!#1]", maxCuts=4, resultsAsMols=True, maxCutBonds=30)
     96 frags += rdMMPA.FragmentMol(mol, pattern="[!#1]!@!=!#[!#1]", maxCuts=3, resultsAsMols=True, maxCutBonds=30)
     97 # hydrogen atoms

RuntimeError: Pre-condition Violation
	RingInfo not initialized
	Violation occurred on line 83 in file Code/GraphMol/RingInfo.cpp
	Failed Expression: df_init
	RDKIT: 2021.09.5
	BOOST: 1_67

This expression works fine when I use a molecule without a dummy atom, for example the core molecule defined above. Any help with this would be greatly appreciated.

Thank you!

Opening Pre-Compiled Databases

Perhaps a novice question, but opening the pre-defined fragment databases given here: http://www.qsar4u.com/pages/crem.php.

using gzip, etc. gives an unexpected end of file error. Any advice on properly unarchiving these?

[bug] package is incompatible with python 3.11+

per the python 3.11 documentation:

The population parameter of random.sample() must be a sequence, and automatic conversion of sets to lists is no longer supported. Also, if the sample size is larger than the population size, a ValueError is raised. (Contributed by Raymond Hettinger in bpo-40465.)

This causes crem to crash when running, as it attempts to use random.sample on a dict.keys object, which apparently is now a non-sequence type (and no longer gets automatically converted).

This is fixed by first converting to a list (or other sequence datatype) prior to calling random.sample.

AtomKekulizeException: non-ring atom 2 marked aromatic

Hello, I am getting a strange error with this molecule.

mol = Chem.MolFromSmiles('C[C@H](C(:[O]):[O-])[C@H](c1ccc2c(c1)O[C@@H](C1CCN([C@@H](C)c3cc(F)ccc3OC(F)(F)F)CC1)CC2)C1CC1')
if mol: print(True)
a =list( mutate_mol(mol, db_name='replacements02_sc2.5.db', radius=3, min_size=0, max_size=2, min_inc=-5, max_inc=6))

The molecule looks like a correct RDKit mol. But i got the following error:

[16:29:26] non-ring atom 2 marked aromatic
---------------------------------------------------------------------------
AtomKekulizeException                     Traceback (most recent call last)
Input In [14], in <cell line: 1>()
----> 1 a =list( mutate_mol(mol, db_name='/home/ale/GITLAB/bi_crem_database/replacements02_sc2.5.db', radius=3, min_size=0, max_size=2, min_inc=-5, max_inc=6))

File ~/anaconda3/envs/moldrug/lib/python3.9/site-packages/crem/crem.py:487, in mutate_mol(mol, db_name, radius, min_size, max_size, min_rel_size, max_rel_size, min_inc, max_inc, max_replacements, replace_cycles, replace_ids, protected_ids, symmetry_fixes, min_freq, return_rxn, return_rxn_freq, return_mol, ncores, **kwargs)
    483 protected_ids = sorted(protected_ids)
    485 if ncores == 1:
--> 487     for frag_sma, core_sma, freq, ids in __gen_replacements(mol1=mol, mol2=None, db_name=db_name, radius=radius,
    488                                                             min_size=min_size, max_size=max_size,
    489                                                             min_rel_size=min_rel_size, max_rel_size=max_rel_size,
    490                                                             min_inc=min_inc, max_inc=max_inc,
    491                                                             max_replacements=max_replacements,
    492                                                             replace_cycles=replace_cycles,
    493                                                             protected_ids_1=protected_ids, protected_ids_2=None,
    494                                                             min_freq=min_freq, symmetry_fixes=symmetry_fixes,
    495                                                             **kwargs):
    496         for smi, m, rxn in __frag_replace(mol, None, frag_sma, core_sma, radius, ids, None):
    497             if max_replacements is None or len(products) < (max_replacements + 1):  # +1 because we added source mol to output smiles

File ~/anaconda3/envs/moldrug/lib/python3.9/site-packages/crem/crem.py:314, in __gen_replacements(mol1, mol2, db_name, radius, dist, min_size, max_size, min_rel_size, max_rel_size, min_inc, max_inc, max_replacements, replace_cycles, protected_ids_1, protected_ids_2, min_freq, symmetry_fixes, **kwargs)
    312 else:
    313     mol = mol1
--> 314     f = __fragment_mol(mol, radius, protected_ids=protected_ids_1, symmetry_fixes=symmetry_fixes)
    316 if f:
    317     mol_hac = mol.GetNumHeavyAtoms()

File ~/anaconda3/envs/moldrug/lib/python3.9/site-packages/crem/crem.py:113, in __fragment_mol(mol, radius, return_ids, keep_stereo, protected_ids, symmetry_fixes)
    110             output.add((env, frag, ids_0))
    111     else:  # multiple cuts
    112         # there are no checks for H needed because H can be present only in single cuts
--> 113         env, frag = get_canon_context_core(chains, core, radius, keep_stereo)
    114         output.add((env, frag, get_atom_prop(core) if return_ids else tuple()))
    116 if symmetry_fixes:

File ~/anaconda3/envs/moldrug/lib/python3.9/site-packages/crem/mol_context.py:302, in get_canon_context_core(context, core, radius, keep_stereo)
    299 def get_canon_context_core(context, core, radius, keep_stereo=False):
    300     # context and core are Mols or SMILES
    301     # returns SMILES by default
--> 302     res = get_std_context_core_permutations(context, core, radius, keep_stereo)
    303     if res:
    304         env, cores = res

File ~/anaconda3/envs/moldrug/lib/python3.9/site-packages/crem/mol_context.py:247, in get_std_context_core_permutations(context, core, radius, keep_stereo)
    245 # remove Hs from context and core
    246 if context:  # context cannot be H (no check needed), if so the user will obtain meaningless output
--> 247     context = Chem.RemoveHs(context)
    248 if core and Chem.MolToSmiles(core) != '[H][*:1]':
    249     core = Chem.RemoveHs(core)

AtomKekulizeException: non-ring atom 2 marked aromatic

Any idea what is going on? A lot of thanks in advance!

Results different between online app and local version

If I draw the following molecule in your web app:

And try to use the "Scaffold" (e.g. grow) functionality at the selected index, I get a number of results (maybe a few hundred/thousand) at radius = 1.

However, if I iterate over the atoms in the molecule in a local version, crem reports it cant generate any results at the same position:

smiles = 'c1cc(OC)ccc1C'
base_mol = Chem.MolFromSmiles(smiles)
atoms = len(list(base_mol.GetAtoms()))

for atom in range(0, atoms):
    gen_results = list(grow_mol(base_mol, return_mol=True, max_replacements=20000,
                                db_name='replacements02_sc2.5.db', radius=1, replace_ids=[atom], max_atoms=15, ncores=4))
    num_results = len(gen_results)
    print(atom, num_results)

The output of this is

[10:36:04] WARNING: not removing hydrogen atom with dummy atom neighbors
0 19999
[10:36:06] WARNING: not removing hydrogen atom with dummy atom neighbors
1 20000
2 0
3 0
[10:36:10] WARNING: not removing hydrogen atom with dummy atom neighbors
4 19999
5 0
6 0
7 0
[10:36:15] WARNING: not removing hydrogen atom with dummy atom neighbors
8 19995
[10:36:18] WARNING: not removing hydrogen atom with dummy atom neighbors
9 19999
[10:36:21] WARNING: not removing hydrogen atom with dummy atom neighbors
10 19998
[10:36:25] WARNING: not removing hydrogen atom with dummy atom neighbors
11 19999
12 0
13 0
14 0
15 0
[10:36:29] WARNING: not removing hydrogen atom with dummy atom neighbors
16 19994
17 0
18 0

And just as a sanity check, this is the numbering of the molecule:

I'm not sure why index 1 is able to produce results, but index 5 doesn't (even though it should by symmetry), and the same position (index 5) works in the web app. Do you know why these would produce different results?

fragmentation py issue: Can't generate mol for: SMILES

Dear Pavel,

I'm so grateful for your powerful work. when I tried to generate the a custom fragment database as you described.
but after runing the py file fragmentation -i input.smi -o frags.txt -c 32 -v, the error happened and indicated that

<Can't generate mol for: SMILES
import sys; print('Python %s on %s' % (sys.version, sys.platform))>

could you please give some suggestions how to fix it and what errors?
neigher did I use the wrong format input.smi file? could you please show me some examples how to use the right smi file?
many many thanks,

best,
Sh-Y

Issue in

Fragmentation of certain molecules fails, it is because ',' as a delimiter doesn't always work . A fix might be to use '\t' as a delimiter as it doesn't appear in SMARTS patterns.

Error in process_line() of frag_to_env_mp.py:

original line CC[C@H]1C[C@@h]1C(=O)NC(C)C,|&1:2,4|,C1C@H[C@H]1[:2],CC(C)NC(=O)[:1].CC[:2]
too many values to unpack (expected 4) ['CC[C@H]1C[C@@h]1C(=O)NC(C)C', '|&1:2', '4|', 'C1C@H[C@H]1[:2]', 'CC(C)NC(=O)[:1].CC[:2]']
There are 5 unpacked values and not 4 as expected because '|&1:2,4|' has a ',' in the pattern

Annoying [13:51:55] WARNING: not removing hydrogen atom with dummy atom neighbors

Hi, excellent package! I just wondering how to stop the warning message that pops up during link_mol and mutate_mol. I tried some

import warnings
warnings.filterwarnings("ignore", message='not removing hydrogen atom with dummy atom neighbors')

without succes.

Problem with guacamol_test.py

Hi all
I was trying the script with guacamol test and got the warning for GPU from tensorflow (I don't have GPU set up), but what I had hard time to understand is this:

(my-rdkit-env) J:\esp\Personal\AndreaZaliani\pubchem_yojana\crem>python crem/guacamol_crem_test.py
2022-07-04 11:22:30.488812: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'cudart64_110.dll'; dlerror: cudart64_110.dll not found
2022-07-04 11:22:30.492783: I tensorflow/stream_executor/cuda/cudart_stub.cc:29] Ignore above cudart dlerror if you do not have a GPU set up on your machine.
Traceback (most recent call last):
File "J:\esp\Personal\AndreaZaliani\pubchem_yojana\crem\crem\guacamol_crem_test.py", line 30, in
from crem import mutate_mol2
File "J:\esp\Personal\AndreaZaliani\pubchem_yojana\crem\crem\crem.py", line 11, in
from crem.mol_context import get_canon_context_core, combine_core_env_to_rxn_smarts
ModuleNotFoundError: No module named 'crem.mol_context'; 'crem' is not a package

Any hint for me?

Best
andrea

error of enumerate_compounds function

Hi, Crem,

I tried the new function enumerate_compounds in utils module.

but there is an error, when I use the method to repeat the example you provied.
73dd185

TerminatedWorkerError: A worker process managed by the executor was unexpectedly terminated. This could be caused by a segmentation fault while calling the function or by an excessive memory usage causing the Operating System to kill the worker.

my PC memory is 32G, so I didnot why this error occurs
could you please provide some suggesions to fix it? many thanks,

Shengyang

auxiliary function get_replacements

how can we implement the function in the new version 0.2.11, auxiliary function get_replacements?

could you please provide some examples how to run the auxiliary function get_replacements? many thanks,

best,

Tutorials

Hi! I was wondering whether you might have some more jupyter tutorials that aren't uploaded on the current GitHub repository? Thanks for building this software, it's super cool!

fragmentation problem

Hi,

I am trying to build my own db according to the instructions, but when I try mols = list(mutate_mol(m, db_name='test.db', max_size=3)) in the tutorial with my own db, the result list of mols is always empty, but when I try with the preproduced db replacements02_sa2.db there is no problem; I used the CHEMBL231.smi file in the example folder, and follow the instructions

fragmentation -i CHEMBL231.smi -o frags.txt -c 32 -v
frag_to_env -i frags.txt -o r3.txt -r 3 -c 32 -v
sort r3.txt | uniq -c > r3_c.txt
env_to_db -i r3_c.txt -o tert.db -r 3 -c -v

got the result

 root@872fabd400c3:/home/crem# python test.py
[]

is there something iI missed?

how to employ the grow_mol2/mut_mol2 function

hi, DrrDom,

could you please provide some example for using grow_mol2 and mut_mol2 functions?
In particular, how to set the parameters, such as replace_ids ?
Many many thanks

enumerate_compounds: filtering the generated molecules to meet my criteria

Hi,

I am using the enumerate_compounds function to generate new molecules, with 3 positions for potential replacement. Hence, I set the n_iterations = 3. As far as I'm aware, this does a stepwise growing of my base molecule mol.

However, I am only interested in the generated molecules where at least 1 specified position is grown. Is there any way to filter out molecules which don't meet this criteria, or how do you suggest I go about this problem?

Thank you so much!

drrdom / crem Goto Github PK

crem's People

Contributors

Stargazers

Watchers

Forkers

crem's Issues

Recommend Projects

Recommend Topics

Recommend Org