Custom QSAR model integration with REINVENT (Staged Learning) not achieving the goal

[filvaleriia]

Dear REINVENT Team,

I am experiencing difficulties integrating my custom QSAR model with REINVENT in staged learning mode. Despite configuring the scoring function and external process, the generated molecules are not achieving the intended optimization goal of predicted pChEMBL values ≥ 7.

Goal:
To create a model that will generate molecules with predicted pChEMBL values are ≥ 7 using my custom QSAR model.

Details:

• Output format for external QSAR scoring function: {"version": 1, "payload": {"predictions": np.reshape(scores, len(smilies)).tolist()}}
• Issue: Only approximately 8% of the generated molecules comply with the condition of pChEMBL ≥ 7, which is significantly lower than expected.

Question:
Could additional adjustments to staged learning configurations help improve alignment with the optimization goal for pChEMBL values?

Below is the code I used for the staged learning setup:

run_type = "staged_learning"
use_cuda = "cuda:0" 
json_out_config = "_staged_learning.json"
output_file = "cOS_REINVENT_staged_learning.csv" 

[parameters]
model_file = "reinvent.prior"

use_checkpoint = false
summary_csv_prefix = "RL"
agent_file = "TL_reinvent.model" #my own model after transfer_learning
prior_file = "reinvent.prior"

batch_size = 512
unique_molecules = true  
randomize_smiles = true  

[learning_strategy]
type = "dap"
sigma = 128
rate = 0.0001

[[stage]]
chkpt_file = "StL.chkpt"

max_score = 1.0
max_steps = 300


[stage.scoring]
type = "geometric_mean"

[[stage.scoring.component]]
[stage.scoring.component.ExternalProcess]

[[stage.scoring.component.ExternalProcess.endpoint]]
name = "QSAR"

min_value = 7
max_value = 10

params.executable = "/usr/bin/conda"
params.args = "run -n env bash -c 'python prediction.py'"

Thank you for your assistance, and please let me know if additional details are needed.

Best regards,
Valeriia Fil

[halx]

Hi,

welcome to the community and many thanks for your interest in REINVENT!

What I note is that min_value and max_value are not valid parameters. In newer versions of REINVENT that would actually be an error. What you probably meant to do here is to provide a transformation function, e.g.

[[stage.scoring.component.ExternalProcess.endpoint]]
name = "QSAR"

params.executable = "/usr/bin/conda"
params.args = "run -n env bash -c 'python prediction.py'"

# Adjust parameters as needed
transform.type = "sigmoid"
transform.high = 10
transform.low = 7
transform.k = 0.4

A step function would probably be too harsh.

Cheers,
Hannes.

[filvaleriia]

Hi,

Thank you very much for your response and for pointing this out! I see now that min_value and max_value are not valid parameters, and it is necessary to use a transformation function instead.

I was wondering if there is a way to guide the model more effectively to generate compounds with a predicted pChEMBL value > 7 (QSAR raw score) or, alternatively, with a transformed score closer to 1? Currently, the majority of generated molecules have predicted values below 7(QSAR (raw)), which is not aligned with the optimization goal.

Looking forward to your advice!

Best regards,
Valeriia Fil

[halx]

Only the transformed score is used in the RL algorithm and needs to be between 0 and 1. You need to make sure to set the transform correctly such that higher pChEMBL values are close to 1 and lower ones closer to 0. That works usually very effectively.

[filvaleriia]

Hello,

Thank you very much for your recommendation regarding the transformation function. I implemented it as suggested, but most of the generated molecules still have a QSAR score close to 0. Consequently, I decided to modify my QSAR model to transform the score internally to a range between 0 and 1 and then used this adjusted model in Staged Learning.

However, despite these adjustments, only about 10% of the molecules achieve a score above 0.6 (the activity threshold). When reviewing my training results, I noticed that the average total score did not exceed 0.48. I have tried modifying various parameters, but these changes have not significantly improved the outcomes.

In the screenshots below, I have included relevant graphs from TensorBoard and the code. Could you kindly provide some guidance or suggestions on how I might adjust my training approach to generate molecules with QSAR scores closer to 1?

Thank you in advance for your time and assistance.

run_type = "staged_learning"
device = "cuda:0"  
tb_logdir = "tb_logs"  
json_out_config = "_staged_learning.json"  

[parameters]
summary_csv_prefix = "RL_dis_0" 
use_checkpoint = false 
purge_memories = false  

## Reinvent
prior_file = "../../priors/reinvent.prior"
agent_file = "TL_reinvent_GR_dis_0.model"

batch_size = 512  

unique_sequences = true  
randomize_smiles = true  
[learning_strategy]

type = "dap"     
sigma = 128       
rate = 0.0001  

[[stage]]

chkpt_file = 'StL_dis_0.chkpt'

termination = "simple"  
max_score = 1.0  
min_steps = 50  
max_steps = 250


[stage.scoring]
type = "geometric_mean" 

[[stage.scoring.component]]
[stage.scoring.component.SAScore]

[[stage.scoring.component.SAScore.endpoint]]
name = "SAScore"
weight = 0.1
transform.type = "sigmoid"
transform.high = 0.1
transform.low = 6
transform.k = 0.5

[[stage.scoring.component]]
[stage.scoring.component.ExternalProcess]

[[stage.scoring.component.ExternalProcess.endpoint]]
name = "QSAR"
weight = 0.9

params.executable = "/usr/bin/conda"
params.args = "run -n env bash -c 'python prediction.py'"

[halx]

Principally I see that there is some optimization going on. But I think it will be difficult to debug this from remote? Would you be able to share your original model?

Also, have you made sure to standardize the SMILES before model building and that the training data set is in-domain with REINVENT? But since you mention pChEMBL, I suppose your model was built from ChEMBL data?

[filvaleriia]

Hi,
Thank you so much for your help! I truly appreciate it. As discussed, I am happy to share all my data, scripts, environment setup, and steps for the project. Below, I’ve outlined the workflow and included relevant code snippets and files:

1. Trasnfer learning: In the first step, I applied transfer learning using my data(cIS_GR_dis_0.smi, cIS_GR_dis_0_valid.smi). I used active compounds targeting the Glucocorticoid receptor. Below is the code snippet for the transfer learning process:

# transfer learning
run_type = "transfer_learning"
device = "cuda:0"
tb_logdir = "tb_TL"
json_out_config = "TL_dis_0.json"
[parameters]
num_epochs = 200 # number of steps to run
save_every_n_epochs = 5 # save checkpoint model file very N steps
batch_size = 100
sample_batch_size = 200 # number of sampled molecules chosen for similarity
## Reinvent
input_model_file = "reinvent.prior"
smiles_file = "cIS_GR_dis_0.smi" # read 1st column
output_model_file = "TL_reinvent_GR_dis_0.model"
validation_smiles_file = "cIS_GR_dis_0_valid.smi"
standardize_smiles = true
randomize_smiles = true
randomize_all_smiles = false
internal_diversity = true

2. Staged learning: After transfer learning, I used the output model (TL_reinvent_GR_dis_0.model) to perform staged learning. The goal was to generate compounds with a SA score < 6 and a pChEMBL value > 7 (or I used a clipped score after prediction to normalize the values between 0 and 1; in this case, I aimed for scores greater than 0.6.). I used a QSAR model from a file, and I’ve attached the necessary files (cIS_Glucocorticoid_receptor_with_p_chembl_dis_0) along with the environment setup (requirements.txt or environment.yml).

Here is the staged learning configuration:

run_type = "staged_learning"
device = "cuda:0"  # set torch device e.g. "cpu"
tb_logdir = "tb_logs"  # name of the TensorBoard logging directory
json_out_config = "_staged_learning.json"  # write this TOML to JSON

[parameters]
summary_csv_prefix = "RL_dis_0"  # prefix for the CSV file
use_checkpoint = false 
purge_memories = false  

## Reinvent
prior_file = "../../priors/reinvent.prior"
agent_file = "TL_reinvent_GR_dis_0.model"

batch_size = 512          # network

unique_sequences = true  
randomize_smiles = true  

[learning_strategy]

type = "dap"      # dap: only one supported
sigma = 128       # sigma of the RL reward function
rate = 0.0001     # for torch.optim

[[stage]]

chkpt_file = 'StL_dis_0.chkpt'

termination = "simple"  # termination criterion fot this stage
max_score = 1.0  # terminate if this total score is exceeded
min_steps = 50  # run for at least this number of steps
max_steps = 250  # terminate entire run when exceeded

[stage.scoring]
type = "geometric_mean"  # aggregation function

[[stage.scoring.component]]
[stage.scoring.component.SAScore]

[[stage.scoring.component.SAScore.endpoint]]
name = "SAScore"
weight = 0.1
transform.type = "sigmoid"
transform.high = 0.1
transform.low = 6
transform.k = 0.5

[[stage.scoring.component]]
[stage.scoring.component.ExternalProcess]

[[stage.scoring.component.ExternalProcess.endpoint]]
name = "QSAR"
weight = 0.9

params.executable = "/home/filv/anaconda3/bin/conda"
params.args = "run --no-capture-output -n env bash -c 'python pchembl_prediction_clipped_score.py'"

3. Sampling: After staged learning, I used the checkpoint to generate molecules. Below is the sampling configuration:

run_type = "sampling"
device = "cuda:0"  # run on the GPU if true, on the CPU if false
json_out_config = "_sampling.json"  # write this TOML to JSON


[parameters]
## Reinvent: de novo sampling
model_file = "StL_dis_0.chkpt"

output_file = 'cOS_REINVENT_dis_0.csv' 
[cIS_Glucocorticoid_receptor_with_p_chembl_dis_0.json](https://github.com/user-attachments/files/18394425/cIS_Glucocorticoid_receptor_with_p_chembl_dis_0.json)
[cIS_Glucocorticoid_receptor_with_p_chembl_dis_0_meta.json](https://github.com/user-attachments/files/18394427/cIS_Glucocorticoid_receptor_with_p_chembl_dis_0_meta.json)
[requirements.txt](https://github.com/user-attachments/files/18394433/requirements.txt)


num_smiles = 1000000
unique_molecules = true  
randomize_smiles = true

Files
I have uploaded all the necessary files to Google Drive. Here’s the link: https://drive.google.com/drive/folders/1zybC_9RZi16Nm1gBH-L9gqW8zIFTn5fe?usp=sharing

Thank you once again for your assistance! Please let me know if you need any additional information or clarification.

Best regards,
Valeriia Fil

[halx]

Hi,

I have pasted am example for staged learning below. Principally, you really need to understand the transformation function and fine-tune accordingly. Sigmoid transforms can have a significant impact on score evaluation (and thus optimization) and may even counter the set weighting factor. Your predictive model for example causes the generated compounds to be in the pChEMBL range of 5.5 to 7.5. Maybe you can push this a bit higher by playing with the transform parameters.

It may also be useful to use a diversity filter to promote structural diversity but more importantly, inception can speed up convergence tremendously. You may also want to consider using a "drug-like" score like QED to at least have some control over the quality of the molecules. SAScore really measures more the complexity of a molecule rather than synthetic tractability. It's a toy score really.

A few other points. TL can make some sense to focus a prior towards known structures before RL but you need to make sure not to overfit. You have a validation set (interestingely quite a bit larger than the actual training set) and if you look at the tensorboard outlook you would see that there is a minimum at around 25 epochs. While I do not principally advocate for using this exact point it is certainly much more sensible than using the final model after 200 steps. The log files even tells you the location of the minimum.

Staged learning (RL) is not typically used to create a model but rather the model is a by-product of the optimization process. "Good" molecules may be generated at any step and you really want to filter (from the CSV file) out those that you think are "good" after the run (you could think about a Pareto approach or any other suitable multi-criteria decision analysis). Post-hoc sampling can't guarantee that the sampling molecules follow your criteria as you do not have the scoring components applied as restraints anymore i.e. if you computed the score on those molecules, many may be uninteresting. Ideally, and often the case, RL will optimize towards a plateau, at which point you are basically sampling under restraints anyway.

There is really no need to run your script through bash. This is just inefficient as it create an unnecessart additional process. You are also recreating the model from JSON on every invocation but the only speed-up here could be gained by transforming that model to a binary format.

Hope that helps,
Hannes.

run_type = "staged_learning"
device = "cuda:0"
tb_logdir = "tb_stage1"
json_out_config = "_stage1.json"


[parameters]

summary_csv_prefix = "stage1"
use_checkpoint = false 
purge_memories = false  

prior_file = ".reinvent"
agent_file = "TL_reinvent_GR_dis_0.model.26.chkpt"

batch_size = 100

unique_sequences = true  
randomize_smiles = true  


[learning_strategy]

type = "dap"
sigma = 128
rate = 0.0001


[inception]

memory_size = 50
sample_size = 10


[diversity_filter]

type = "IdenticalMurckoScaffold"
bucket_size = 10
minscore = 0.7
minsimilarity = 0.5


[[stage]]

chkpt_file = 'stage1.chkpt'

termination = "simple"
max_score = 1.0
max_steps = 500


[stage.scoring]
type = "geometric_mean"  # aggregation function


[[stage.scoring.component]]
[stage.scoring.component.SAScore]

[[stage.scoring.component.SAScore.endpoint]]
name = "SAScore"
weight = 0.1
transform.type = "reverse_sigmoid"
transform.high = 10.0
transform.low = 1.0
transform.k = 0.5


[[stage.scoring.component]]
[stage.scoring.component.QED]

[[stage.scoring.component.QED.endpoint]]
name = "QED"
weight = 0.2


[[stage.scoring.component]]
[stage.scoring.component.ExternalProcess]

[[stage.scoring.component.ExternalProcess.endpoint]]
name = "QSAR"
weight = 0.7

params.executable = "/home/kmll395/miniconda3/bin/conda"
params.args = "run -n qsprpred --no-capture-output /home/kmll395/projects/reinvent/public/filvaleriia/pchembl_prediction.py"  # no bash please!

transform.type = "sigmoid"
transform.high = 7.5
transform.low = 5.5
transform.k = 0.5

[halx]

Oh, and and a basic alerts filter to avoid some of the more offensive compounds (mid-sized rings, typically).

[[component]]
[component.custom_alerts]

[[component.custom_alerts.endpoint]]
name = "Alerts"

params.smarts = [
    "[*;r8]",
    "[*;r9]",
    "[*;r10]",
    "[*;r11]",
    "[*;r12]",
    "[*;r13]",
    "[*;r14]",
    "[*;r15]",
    "[*;r16]",
    "[*;r17]",
    "[#8][#8]",
    "[#6;+]",
    "[#16][#16]",
    "[#7;!n][S;!$(S(=O)=O)]",
    "[#7;!n][#7;!n]",
    "C#C",
    "C(=[O,S])[O,S]",
    "[#7;!n][C;!$(C(=[O,N])[N,O])][#16;!s]",
    "[#7;!n][C;!$(C(=[O,N])[N,O])][#7;!n]",
    "[#7;!n][C;!$(C(=[O,N])[N,O])][#8;!o]",
    "[#8;!o][C;!$(C(=[O,N])[N,O])][#16;!s]",
    "[#8;!o][C;!$(C(=[O,N])[N,O])][#8;!o]",
    "[#16;!s][C;!$(C(=[O,N])[N,O])][#16;!s]"
]