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: 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 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..3a3bfd2 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, + 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, + 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) @@ -176,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 15bc960..e06b12f 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 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 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: @@ -265,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( @@ -273,11 +298,11 @@ 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 - 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 @@ -312,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 @@ -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'])