From c679055a6a579a854c79f72775451616e71bcc4a Mon Sep 17 00:00:00 2001 From: MariBerry Date: Tue, 3 Mar 2026 09:17:24 +0100 Subject: [PATCH 1/5] added bounding_box column to output database --- optimizer.py | 8 ++++-- optimizer_utils.py | 38 ++++++++++++++++++-------- process_predictions.py | 62 ++++++++++++++++++++++++++++++------------ 3 files changed, 77 insertions(+), 31 deletions(-) diff --git a/optimizer.py b/optimizer.py index 177b5fa..0fb569c 100644 --- a/optimizer.py +++ b/optimizer.py @@ -30,10 +30,13 @@ def optimize(settings: Dict, input_config: str, brute_force: bool) -> None: parameters_list_dicts = [settings[parameter] for parameter in settings if "param_" in parameter] + use_spci_models = settings['descriptors_type'] != "MPNN_fingerprint" + # create database settings['output_database'] = optimizer_utils.create_database( settings['working_dir'], - [parameter['name'] for parameter in parameters_list_dicts] + [parameter['name'] for parameter in parameters_list_dicts], + use_spci_models ) shutil.copyfile(input_config, os.path.join(settings['working_dir'], 'config.yaml')) @@ -156,7 +159,8 @@ def optimize(settings: Dict, input_config: str, brute_force: bool) -> None: desirabilities, settings['number_of_selected_compounds'], settings['random_compounds_selection'], - brute_force + brute_force=brute_force, + spci_models=use_spci_models ) # get num of fitted compounds diff --git a/optimizer_utils.py b/optimizer_utils.py index 593e50c..35569a8 100644 --- a/optimizer_utils.py +++ b/optimizer_utils.py @@ -49,28 +49,42 @@ def save_output_poll(in_sdf: str, out_fname: str, output_poll: pandas_table) -> str(output_poll.loc[mol.GetProp('id'), col])) output.write(mol) -def create_database(working_dir: str, parameter_to_optimize: List) -> str: +def create_database(working_dir: str, parameter_to_optimize: List, spci_models: bool = False) -> str: """ Define new database and save it to working_dir. :param working_dir: path to working directory, new db will be stored there :param parameter_to_optimize: list of all parameters + :param spci_models: bool specifying whether spci models are used (and thus bounding_box should be saved) :return: path_to_database """ parameter_to_optimize = ["predicted_{}".format(parameter) for parameter in parameter_to_optimize] path_to_database = os.path.join(working_dir, 'output.db') - table_str = ("""CREATE TABLE optimizer_table( - id TEXT NOT NULL, - smi TEXT NOT NULL UNIQUE, - mol_block TEXT NOT NULL UNIQUE, - protected_ids TEXT NOT NULL, - generation INTEGER NOT NULL, - parent TEXT, - transformation TEXT, - fit INTEGER, """ + - " REAL,".join(parameter_to_optimize) + - " REAL)") + if spci_models: + table_str = ("""CREATE TABLE optimizer_table( + id TEXT NOT NULL, + smi TEXT NOT NULL UNIQUE, + mol_block TEXT NOT NULL UNIQUE, + protected_ids TEXT NOT NULL, + generation INTEGER NOT NULL, + parent TEXT, + transformation TEXT, + fit INTEGER, """ + + " REAL,".join(parameter_to_optimize) + + " REAL, bounding_box INTEGER)") + else: + table_str = ("""CREATE TABLE optimizer_table( + id TEXT NOT NULL, + smi TEXT NOT NULL UNIQUE, + mol_block TEXT NOT NULL UNIQUE, + protected_ids TEXT NOT NULL, + generation INTEGER NOT NULL, + parent TEXT, + transformation TEXT, + fit INTEGER, """ + + " REAL,".join(parameter_to_optimize) + + " REAL)") if os.path.isfile(path_to_database): os.remove(path_to_database) diff --git a/process_predictions.py b/process_predictions.py index 15bc960..1858688 100644 --- a/process_predictions.py +++ b/process_predictions.py @@ -74,7 +74,7 @@ def save_output(input_sdf: str, out_fname: str, output_poll: pandas_table) -> No in_file.close() -def prepare_working_arr(in_pred: List, parameters: List, bounding_box: bool, proba_consensus:bool=True) -> pandas_table: +def prepare_working_arr(in_pred: List, parameters: List, bounding_box: bool, proba_consensus:bool=True, spci_models:bool=False) -> pandas_table: # todo param proba_consensus should go to config, so regression recalculation will be avoided+bettertracking of run """ Reads file with predictions and process it into pandas table @@ -98,17 +98,34 @@ def prepare_working_arr(in_pred: List, parameters: List, bounding_box: bool, pro else: table.drop(table[table.bound_box==0].index, inplace=True) - table.drop('bound_box', axis=1, inplace=True) if proba_consensus: - table.drop(table.columns[-1], axis=1, inplace=True)# drop consensus - table['consensus'] = table.loc[:,['Compounds' not in i for i in table.columns]].mean(axis=1) # get new consensus # TODO: PP, why "consensus" was used? why not to remove the first column by index? - - cols_to_drop = list(range(1,table.shape[1]-1)) # drop all but (new) consensus - table.drop(table.columns[cols_to_drop], axis=1, inplace=True) - table.rename(columns={table.columns[0]: 'id', 'consensus': parameter}, inplace=True) + if 'consensus' in table.columns: + table.drop('consensus', axis=1, inplace=True) + model_cols = [c for c in table.columns if 'Compounds' not in str(c) and c != 'bound_box'] + table['consensus'] = table[model_cols].mean(axis=1) + + cols_to_keep = [table.columns[0], 'consensus'] + if spci_models and 'bound_box' in table.columns: + cols_to_keep.append('bound_box') + + cols_to_drop = [c for c in table.columns if c not in cols_to_keep] + table.drop(cols_to_drop, axis=1, inplace=True) + + rename_dict = {table.columns[0]: 'id', 'consensus': parameter} + if spci_models and 'bound_box' in table.columns: + rename_dict['bound_box'] = f'bound_box_{parameter}' + + table.rename(columns=rename_dict, inplace=True) table.set_index('id', inplace=True) tables = pd.concat(tables, axis=1, join='inner') + + if spci_models: + bb_cols = [col for col in tables.columns if str(col).startswith('bound_box_')] + if bb_cols: + # Bounding box is identical for all SPCI models, no need to take minimum value + tables['bounding_box'] = tables[bb_cols[0]].astype(int) + tables.drop(bb_cols, axis=1, inplace=True) # if ad if bounding_box: # TODO: PP, this condition is not needed @@ -139,7 +156,7 @@ def compute_distance_from_threshold(x: float, threshold: List) -> float: def update_database(out_database: str, predictions: pandas_table, - output_filtering: pandas_table) -> None: + output_filtering: pandas_table, spci_models: bool = False) -> None: """ Updates predicted values for compounds in database. @@ -148,11 +165,15 @@ def update_database(out_database: str, predictions: pandas_table, :param output_filtering: fitted compounds in table """ - columns = predictions.columns + columns = [col for col in predictions.columns if col != 'bounding_box'] database_columns = ["predicted_{}".format(cols) for cols in columns] - query = "UPDATE optimizer_table SET fit=?, " - query += "=?, ".join(database_columns) + "=? WHERE id=?" + if spci_models and 'bounding_box' in predictions.columns: + query = "UPDATE optimizer_table SET fit=?, " + query += "=?, ".join(database_columns) + "=?, bounding_box=? WHERE id=?" + else: + query = "UPDATE optimizer_table SET fit=?, " + query += "=?, ".join(database_columns) + "=? WHERE id=?" con = lite.connect(out_database) with con: @@ -166,6 +187,10 @@ def update_database(out_database: str, predictions: pandas_table, record.append(0) for col in columns: record.append(row[col]) + + if spci_models and 'bounding_box' in predictions.columns: + record.append(int(row['bounding_box'])) + record.append(index) cursor.execute(query, tuple(record)) @@ -213,7 +238,7 @@ def get_norm_value(x_input: float, function: List) -> float: def main(in_sdf, in_pred, out_database, out_fname, parameters, optimization_method, thresholds, ad, desirabilities=None, - n_compounds=0, random_compounds=0, additive_agg = True, brute_force=False): + n_compounds=0, random_compounds=0, additive_agg = True, brute_force=False, spci_models=False): """ Logic of algorithm: 1. if brute force - save all compounds to out_fname. @@ -241,11 +266,11 @@ def main(in_sdf, in_pred, out_database, out_fname, parameters, print('Processing predictions ...') # process all predictions - predictions = prepare_working_arr(in_pred, parameters, ad) + predictions = prepare_working_arr(in_pred, parameters, ad, spci_models=spci_models) if ad and predictions is None: print('Compounds are not in ad. Calculating outside ad!') - predictions = prepare_working_arr(in_pred, parameters, False) + predictions = prepare_working_arr(in_pred, parameters, False, spci_models=spci_models) # check if brute_force is selected, then no selections if brute_force: @@ -273,7 +298,7 @@ def main(in_sdf, in_pred, out_database, out_fname, parameters, os.path.join(os.path.dirname(in_sdf), 'output_match.sdf'), output_filtering ) - update_database(out_database, prepare_working_arr(in_pred, parameters, False), output_filtering) + update_database(out_database, prepare_working_arr(in_pred, parameters, False, spci_models=spci_models), output_filtering, spci_models=spci_models) # calc random compounds needed number n_random = math.floor(n_compounds * random_compounds) # will be used later @@ -370,9 +395,12 @@ def main(in_sdf, in_pred, out_database, out_fname, parameters, parser.add_argument('-bf', '--brute_force', action='store_true', default=False, help='use all compounds, no selections') + parser.add_argument('-spci', '--spci_models', action='store_true', default=False, + help='if spci models are used') + args = vars(parser.parse_args()) main(args['in_sdf'], args['in_pred'], args['output_database'], args['out'], args['parameters'], args['methods'], args['thresholds'], args['ad'], args['desirabilities'], args['n_compounds'], args['random_compounds'], - args['brute_force']) + additive_agg=True, brute_force=args['brute_force'], spci_models=args['spci_models']) From 586c58ef184002d7d6b1dcef67b526a0d7049691 Mon Sep 17 00:00:00 2001 From: MariBerry Date: Tue, 31 Mar 2026 16:55:22 +0200 Subject: [PATCH 2/5] bugfix --- process_predictions.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/process_predictions.py b/process_predictions.py index 1858688..9004adb 100644 --- a/process_predictions.py +++ b/process_predictions.py @@ -290,7 +290,7 @@ def main(in_sdf, in_pred, out_database, out_fname, parameters, threshold=threshold) # find compounds which are in threshold - output_filtering = distance_predictions[distance_predictions.apply(lambda x: np.all(x<=0), axis=1)] # all parameteres within thres + output_filtering = distance_predictions[distance_predictions[parameters].apply(lambda x: np.all(x<=0), axis=1)] # all parameteres within thres output_filtering = predictions.loc[output_filtering.index].copy() if output_filtering.shape[0] > 0: save_output( @@ -302,7 +302,7 @@ def main(in_sdf, in_pred, out_database, out_fname, parameters, # calc random compounds needed number n_random = math.floor(n_compounds * random_compounds) # will be used later - ids_not_in_thr = distance_predictions.apply(lambda x: np.any(x > 0), axis=1) # ids of compounds not in threshold + ids_not_in_thr = distance_predictions[parameters].apply(lambda x: np.any(x > 0), axis=1) # ids of compounds not in threshold # print(ids_not_in_thr) # if we have less or equal compounds in input sdf then we specified @@ -337,9 +337,9 @@ def main(in_sdf, in_pred, out_database, out_fname, parameters, apply(get_norm_value, function=function) if additive_agg: # additive - desirability_predictions['desirability'] = desirability_predictions.sum(axis=1)/(len(parameters)) + desirability_predictions['desirability'] = desirability_predictions[parameters].sum(axis=1)/(len(parameters)) else: # multiplicative - desirability_predictions['desirability'] = desirability_predictions.product(axis=1) + desirability_predictions['desirability'] = desirability_predictions[parameters].product(axis=1) desirability_predictions = desirability_predictions.sort_values(by='desirability', ascending=False) selected_compounds_index = predictions.loc[desirability_predictions.head(n_compounds).index].index From f79c3196cc5c8b1fd014aaf5a913f0d45e35a27e Mon Sep 17 00:00:00 2001 From: MariBerry Date: Wed, 8 Apr 2026 12:12:51 +0200 Subject: [PATCH 3/5] update examples --- .../input_config.yaml | 18 +++++++++--------- .../readme.txt | 4 ++-- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/examples/antihistamine_drug_optimization_MPNN/input_config.yaml b/examples/antihistamine_drug_optimization_MPNN/input_config.yaml index 5b28158..b73c302 100644 --- a/examples/antihistamine_drug_optimization_MPNN/input_config.yaml +++ b/examples/antihistamine_drug_optimization_MPNN/input_config.yaml @@ -1,37 +1,37 @@ -working_dir: examples/antihistamine_drug_optimization_SPCI_SCIKIT_LEARN +working_dir: examples/antihistamine_drug_optimization_MPNN number_of_selected_compounds: 10 seed_structure: examples/antihistamine_drug_optimization_SPCI_SCIKIT_LEARN/decloxizine_seed.sdf num_of_generations: 5 num_output_compounds: 10 n_cores: 1 random_compounds_selection: 0 -optimization_method: desirability +optimization_method: pareto smarts_string: '[!#1]!@!=!#[!#1]' max_cuts: 1 protected_ids: 'protected_ids' replacement_database: examples/replacements02_sa2.db -number_of_worst_fragments: 6 +number_of_worst_fragments: 8 random_fragments_selection: 0 max_frag_size: 7 min_inc: -2 max_inc: 2 radius: 2 -descriptors_type: AP +descriptors_type: MPNN_fingerprint bounding_box: True multitask: False param_0: name: 231 - path: examples/models/231_spci_sklearn_atompairs_fp/ - types_of_alg: rf gbm + path: examples/models/231_checkpoints/ + types_of_alg: MPNN type_of_model: reg threshold: more7.2 range: 6.3 desirability: 6:0,10:0.25*x-1.5,100:1 param_1: name: logbb - path: examples/models/logbb_spci_sklearn_atompairs_fp/ - types_of_alg: rf gbm + path: examples/models/logbb_checkpoints/ + types_of_alg: MPNN type_of_model: class threshold: less0.5 range: 1 - desirability: 0:1,1:1-x,100:0 \ No newline at end of file + desirability: 0:1,1:1-x,100:0 diff --git a/examples/antihistamine_drug_optimization_MPNN/readme.txt b/examples/antihistamine_drug_optimization_MPNN/readme.txt index 033b351..f3eca11 100755 --- a/examples/antihistamine_drug_optimization_MPNN/readme.txt +++ b/examples/antihistamine_drug_optimization_MPNN/readme.txt @@ -1,5 +1,5 @@ -Please, download replacement databse and replace path in config file, line 15: +Please, download replacement databse and replace path in config file: replacement_database: examples/replacements02_sa2.db @@ -7,4 +7,4 @@ To run this example, type in command line $ cd crem-ml $ conda activate crem_ml_env -$ python optimizer.py -i examples/antihistamine_drug_optimization_MPNN/input_config.yaml \ No newline at end of file +$ python optimizer.py -i examples/antihistamine_drug_optimization_MPNN/input_config.yaml From 3ae7ffac97d8eb8d19ccad4b1d2753fbf3184792 Mon Sep 17 00:00:00 2001 From: MariBerry Date: Sat, 11 Apr 2026 16:39:41 +0200 Subject: [PATCH 4/5] fix handling of potected_ids when they are not provided or None provided.fixed bounding_box to use all parameters not just one, they all have different AD, choose using logical AND --- optimizer_utils.py | 7 ++++--- process_predictions.py | 4 ++-- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/optimizer_utils.py b/optimizer_utils.py index 35569a8..3a3bfd2 100644 --- a/optimizer_utils.py +++ b/optimizer_utils.py @@ -66,7 +66,7 @@ def create_database(working_dir: str, parameter_to_optimize: List, spci_models: id TEXT NOT NULL, smi TEXT NOT NULL UNIQUE, mol_block TEXT NOT NULL UNIQUE, - protected_ids TEXT NOT NULL, + protected_ids TEXT, generation INTEGER NOT NULL, parent TEXT, transformation TEXT, @@ -78,7 +78,7 @@ def create_database(working_dir: str, parameter_to_optimize: List, spci_models: id TEXT NOT NULL, smi TEXT NOT NULL UNIQUE, mol_block TEXT NOT NULL UNIQUE, - protected_ids TEXT NOT NULL, + protected_ids TEXT, generation INTEGER NOT NULL, parent TEXT, transformation TEXT, @@ -190,7 +190,8 @@ def get_mols(database: str, gen: int = None, fields: list = None) -> list: mol = Chem.MolFromMolBlock(item[0], removeHs=False) if mol: for field, value in zip(fields, item[1:]): - mol.SetProp(field, value) + if value is not None: + mol.SetProp(field, value) mols.append(mol) return mols diff --git a/process_predictions.py b/process_predictions.py index 9004adb..e06b12f 100644 --- a/process_predictions.py +++ b/process_predictions.py @@ -123,8 +123,8 @@ def prepare_working_arr(in_pred: List, parameters: List, bounding_box: bool, pro if spci_models: bb_cols = [col for col in tables.columns if str(col).startswith('bound_box_')] if bb_cols: - # Bounding box is identical for all SPCI models, no need to take minimum value - tables['bounding_box'] = tables[bb_cols[0]].astype(int) + # Bounding box is evaluated as the minimum (AND operation) across all SPCI models + tables['bounding_box'] = tables[bb_cols].min(axis=1).astype(int) tables.drop(bb_cols, axis=1, inplace=True) # if ad From a14bf5e9fb1eb3fd1ec4cde710caa3db1e6cb2a6 Mon Sep 17 00:00:00 2001 From: MariBerry Date: Tue, 12 May 2026 12:57:28 +0200 Subject: [PATCH 5/5] updated README:md --- README.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 3e20c8f..2f3f510 100644 --- a/README.md +++ b/README.md @@ -153,7 +153,7 @@ or pareto (not both). Desirability ranks compounds according to desirability fun lead to these hydrogens being selected as worst fragments, because of 0 contributions. + max_cuts - number of maximum cuts used in fragmentation procedure (see RDKit.Chem.rdMMPA docs). Recommended is 1...3 cuts. + protected_ids - field name in sdf file (with structure(s) to be optimized), which contains 0-based ids of atoms, - that should not be affected by optimization. For instance, this could be an active scaffold which you wish to retain. + that should not be affected by optimization. The format should be field name in SDF-file with the seed molecule(s) containing a comma-separated 1-bassed atom ids. For instance, this could be an active scaffold which you wish to retain. + replacement_database - path to database with interchangeable fragments. NOTE: If invalid database is supplied, no compounds will be generated, and next steps can result in errors (without warning). @@ -226,6 +226,8 @@ TT (topological torsion); or MPNN_fingerprint. Prefix b- means binary fingerprin Building models before running CReM-ML -------------------------------------- +**IMPORTANT:** SPCI models MUST be built using a training set SDF file with explicit hydrogens. This is critical because CReM-ML uses explicit hydrogens during the optimization process. Furthermore, for information regarding which descriptor types are eligible for use with SPCI and how they are encoded, please refer to the command line help message of `spci_descriptors` (or `descriptors.py`) in the SPCI package. + CReM-ML requires ready QSAR models for properties to be optimized. Compatible models are scikit-learn models and Chemprop MPNN models (pytorch-based). Scikit-learn models can be built using SPCI package using SPCI GUI: