added bounding_box column to output database

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:

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
+  desirability: 0:1,1:1-x,100:0

examples/antihistamine_drug_optimization_MPNN/readme.txt

@@ -1,10 +1,10 @@
 
-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
 
 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
+$ python optimizer.py -i examples/antihistamine_drug_optimization_MPNN/input_config.yaml

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

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
 

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,19 +290,19 @@ 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(
                 in_sdf,
                 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'])