GenCReM: molecule optimization, scaffold decoration, de novo generation guided by explainable docking
- Molecule optimization, scaffold decoration, de novo generation tasks
- Based on a genetic algorithm (GA)
- Uses the CReM[1] engine to generate feasible, synthetically accessible molecules
- Makes it easy to specify physicochemical property limits for generated molecules
- Optimizes docking score, number of protein-ligand interactions, and RMSD
- Indirectly controls synthetic accessibility via CReM parameters
- Explores chemical space via an explainable docking workflow that preserves well-fitted 3D fragments while editing suboptimal atoms
Dependencies
- Python (>=3.11)
conda env create --file environment.yaml
conda activate gencrem
pip install git+https://github.com/drrdom/crem.git
pip install git+https://github.com/ci-lab-cz/easydock.git
git clone https://github.com/ci-lab-cz/crem-dock-ga.git
For protonation, install the required packages (MolGPKA, UNI-PKA, etc.) from the easydock installation page.
Note
Fragment databases are required for the CReM generator and can be found here.
python gencrem.py -h
usage: gencrem.py [-h] [-i FILENAME] -o FILENAME [--hostfile FILENAME] [--log FILENAME] [--prefix STRING] [-c NCPU] [--seed SEED] [-d fragments.db] [-r RADIUS]
[--min_freq MIN_FREQ] [--replace_cycles] [--max_size MAX_SIZE] [--max_inc MAX_INC] [--min_inc MIN_INC] [--n_iteration N_ITERATION]
[--population_size POPULATION_SIZE] [--population_replicas {1}] [--mutation_rate MUTATION_RATE] [--elitism] [--minimize_fitness]
[--max_lifetime MAX_LIFETIME] [--mode {parallel,survival}] [--norm_func {softmax,piecewise}] [--stop_if_converged]
[--convergence_rate CONVERGENCE_RATE] [--convergence_method {duplication,improvement}] [--no_duplicates] [--no_contrib] [--rmsd RMSD] [--mw MW]
[--rtb RTB] [--logp LOGP] [--tpsa TPSA] [--plif [PLIF ...]] [--plif_cutoff NUMERIC]
[--plif_map [scaled values for each PLIF similarity value: 0 contacts, 1 contact ... any Nth contacts ...]] [--no_tautomerization] [--no_redock]
[--protonation PROTONATION] [--protein_h protein.pdb] [--program STRING] [--config FILENAME]
Fragment mutation within the binding pocket using Docking contributions.
options:
-h, --help show this help message and exit
-i FILENAME, --input_frags FILENAME
SMILES file with input fragments/molecules or SDF file with 3D coordinates of pre-aligned input fragments (e.g. from PDB complexes). If SDF contains
a <protected_user_ids> field (comma-separated 1-based indices) these atoms will be protected from mutation. This argument can be omitted if an
existing output DB is specified; docking will then continue from the last successful iteration. Optional. (default: None)
-o FILENAME, --output FILENAME
SQLite DB with docking results. If an existing DB was supplied, input fragments (if any) will be ignored and the program will continue docking
from the last successful iteration. (default: None)
--hostfile FILENAME text file with addresses of nodes of a Dask SSH cluster. Most typically, it can be passed as the $PBS_NODEFILE variable from inside a PBS script.
The first line in this file will be the address of the scheduler running on the standard port 8786. If omitted, calculations run on a
single machine as usual. (default: None)
--log FILENAME log file to collect progress and debug messages. If omitted, the log file with the same name as output DB will be created. (default: None)
--prefix STRING prefix to add to all names; useful if multiple runs are analyzed together. (default: None)
-c NCPU, --ncpu NCPU number of CPUs. (default: 24)
--seed SEED Seed for the Random class (default: None)
--population_size POPULATION_SIZE
Size of replica populations (default: 20)
--no_duplicates Do not allow generating duplicate candidates that were already generated (default: False)
--no_contrib Do not use atom contribution (default: False)
CReM parameters:
-d fragments.db, --db fragments.db
SQLite DB with fragment replacements. (default: None)
-r RADIUS, --radius RADIUS
context radius for replacement. (default: 3)
--min_freq MIN_FREQ frequency of occurrence of the fragment in the source database. (default: 0)
--replace_cycles look for a replacement of a fragment containing cycles regardless of fragment size. (default: False)
--max_size MAX_SIZE maximum number of heavy atoms in a fragment to replace. By default, calculated automatically based on docking score. The user preset value will
be ignored for molecules with no replacements found with such parameters. (default: None)
--max_inc MAX_INC maximum change in the number of heavy atoms in replacing fragments relative to the number in the replaced one. By default, calculated
automatically based on MW. The user preset value will be ignored for molecules with no replacements found with such parameters. (default: None)
--min_inc MIN_INC default minimum change in the number of heavy atoms in replacing fragments relative to the number in the replaced one. Negative values mean the
replacing fragments are smaller than the replaced one by the specified number of heavy atoms. By default, calculated automatically based
on MW. The user preset value will be ignored for molecules with no replacements found with such parameters. (default: None)
Genetic parameters:
--n_iteration N_ITERATION
Number of generations to evolve. If you want to continue from the previous last iteration (output DB argument should match the previous run),
set the total number of generations you want to perform. (default: 5)
--mutation_rate MUTATION_RATE
Rate of mutation operation (default: 0.7)
--elitism Preserve the best-fitting molecule for the new replicas populations (default: False)
--minimize_fitness Minimize fitness. Default is to maximize. (default: False)
--max_lifetime MAX_LIFETIME
To prevent population stagnation, restart a molecule to the initial molecule if it stays unchanged for 20 generations in a row. (default: 20)
--mode {parallel,survival}
GenCReM optimization mode (default: parallel)
--norm_func {softmax,piecewise}
GenCReM fitness normalization function (default: softmax)
--stop_if_converged Stop the run when, for the last 5 iterations, the rate of duplicated molecules is >= 0.85 (this value can be changed), meaning the algorithm
has already explored local space and cannot find more new solutions. (default: False)
--convergence_rate CONVERGENCE_RATE
Stop the run if --stop_if_converged is used and the rate of duplicated molecules is >= the --convergence_rate value for the last 5 iterations.
(default: 0.85)
--convergence_method {duplication,improvement}
Stop the run if --stop_if_converged is used. Duplication method: stop if for the last 5 iterations the duplicate rate in the population
is higher than the --convergence_rate argument value. Improvement method: stop if for the last 5 iterations no improvement of docking score was shown.
(default: duplication)
Additional molecule parameters:
--rmsd RMSD maximum allowed RMSD value relative to a compound to pass to the next iteration. By default, no filtering (default: None)
--mw MW maximum ligand molecular weight to pass on the next iteration. (default: 450)
--rtb RTB maximum allowed number of rotatable bonds in a compound. (default: 5)
--logp LOGP maximum allowed logP of a compound. (default: 4)
--tpsa TPSA maximum allowed TPSA of a compound. (default: 120)
--plif [PLIF ...] list of protein-ligand interactions compatible with ProLIF. Dot-separated names of each interaction which should be observed for a ligand.
Derive these names from a reference ligand. Example: ASP115.HBDonor or ARG34.A.Hydrophobic. Choose interaction type from [Hydrophobic,
HBDonor, HBAcceptor, Anionic, Cationic, CationPi, PiCation, FaceToFace, EdgeToFace, MetalAcceptor] (default: None)
--plif_cutoff NUMERIC
cutoff of Tversky similarity, value between 0 and 1. (default: 1)
--plif_map [scaled values for each PLIF similarity value: 0 contacts, 1 contact ... any Nth contacts ...]
Example: 0 0.5 0.9 1 (default: None)
Easydock parameters:
--no_tautomerization disable tautomerization of molecules during protonation (applicable to Chemaxon only). (default: False)
--no_redock Do not run redocking for the input SDF conformation. Use the input conformation. (default: False)
--protonation PROTONATION
Choose a protonation program supported by EasyDock. (default: None)
--protein_h protein.pdb
PDB file with the same protein as for docking, but it should have all hydrogens explicit. Required for PLIF detection. (default: None)
--program STRING name of a docking program supported by EasyDock
--config FILENAME YAML file with parameters used by the docking program (e.g., vina.yml or gnina.yml). vina.yml: protein path to a PDBQT file with a protein;
protein_setup: path to a text file with binding-site coordinates. (default: None)
gencrem.py -i example/2btr_A_ligand_protected_atoms.sdf -o output.db \
-d chembl33_f5.db --protein_h example/2btr_A_H.pdb \
--config example/config_easydock_gnina.yaml \
--population_size 10 --max_size 10 \
--seed 1024 --mutation_rate 1 \
--n_iteration 100 --no_tautomerization \
--program gnina --protonation molgpka \
--radius 3 --no_duplicates \
--mode survival --elitism --norm_func softmax \
--stop_if_converged --max_lifetime 15 \
--convergence_rate 0.9 --convergence_method duplication \
--no_redock
To force GenCReM to preserve important contacts, use additional arguments:
--plif leu83.a.hbdonor leu83.a.hbacceptor --plif_cutoff 0.5 --plif_map 0 0.5 0.9
- Input SDF should not contain hydrogens, or hydrogens should be listed after heavy atoms; otherwise protected IDs can be broken and will not be correctly transferred to child molecules.
- List atom indices to protect as a
protected_user_idsproperty in the SDF file for each molecule (seeexample/2btr_A_ligand_protected_atoms.sdf). - Indices of preserved atoms should start from 1.
- Example SDF molecule property:
<protected_user_ids>
3,4,5,6,7,8,19
Early stop when no new molecules are found for the last N iterations (5 by default), which means the algorithm has already explored local space and cannot find more new solutions
--stop_if_converged - to use preliminary stop mode
--convergence_rate 0.9 - if 90% of new molecules are seen already in the previous iteration
--convergence_method duplication/improvement (default: duplication)
Methods:
duplication - stop if for the last --n_iter_converged iterations the duplicate rate in the population is higher than the --convergence_rate argument value.
improvement - stop if for the last --n_iter_converged iterations no improvement of docking score was shown.
--n_iter_converged 5 - stop if for the last 5 iterations convergence is reached
Reset a molecule to a random initial seed molecule if it has not improved for --max_lifetime iterations in a row.
- Polishchuk, P., CReM: Chemically reasonable mutations framework for structure generation, Journal of Cheminformatics, 2020, 12, 1, 1-18.