Category: Uncategorized

Mol2Mol in REINVENT
I recently attended UK QSAR event and Keishi Kohara from AstraZeneca gave a great talk on REINVENT4, using Mol2Mol transformer. I was vaguely aware of the option but had never used it. From Keishi’s talk it looked a lot simpler than the other “whack-a-molecule” options of iteratively editing input TOML files. Here is the relevant description in the paper.

Molecular optimization [22, 52]. A molecule is supplied to the generator as restraint. The generator will find a second molecule within a defined similarity. Depending on the similarity radius the molecule will be relatively similar to the supplied molecule but, importantly, the scaffold can change within the limits of the given similarity. (Mol2Mol)

Basically you just pick a pre-trained prior, a list of project molecules along with a few details and sample the model.

Install

Install is really easy, switch runtime to gpu/tpu, check what version of cuda is running (in this case 12.5). Set up conda for colab.
```
!nvcc --version
!pip install -q condacolab
import condacolab
condacolab.install_miniforge()
```
```
!git clone https://github.com/MolecularAI/REINVENT4.git
%cd REINVENT4
!python install.py cu125
```
That’s it, I wish installing other python packages was so easy! Next step is writing a toml file. We’ll create analogues around OSIMERTINIB where the smiles is saved in “osimertinib.smi”.
```
# REINVENT4 TOML input example for sampling
#

run_type = "sampling"
device = "cuda:0"  # set torch device e.g. "cpu"
json_out_config = "_sampling.json"  # write this TOML to JSON

[parameters]

## Mol2Mol: find molecules similar to the provided molecules
model_file = "priors/mol2mol_medium_similarity.prior"
smiles_file = "osimertinib.smi"  # 1 compound per line
sample_strategy = "beamsearch"  # multinomial or beamsearch (deterministic)
temperature = 1.0 # temperature in multinomial sampling
#tb_logdir = "tb_logs"  # name of the TensorBoard logging directory

output_file = 'sampling.csv'  # sampled SMILES and NLL in CSV format

num_smiles = 100000  # number of SMILES to be sampled, 1 per input SMILES
unique_molecules = true  # if true remove all duplicatesd canonicalize smiles
randomize_smiles = true # if true shuffle atoms in SMILES randomly
```
Running the code is very simple
```
!reinvent -l sampling.log sampling.toml
```
7th May 2025

BDE with orca compound script

Most of the time for Orca QM calculations I’ll use Weasel FACCT’s rather fantastic workflow tool but occasionally there are workflows which either don’t exist or don’t quite do exactly what I want. After attending the Virtual Winter School on Computational Chemistry and getting lessons on Dimitrios Liakos on compound scripts I’ve been mean to put the theory into practice. The talks on line here , if your new to Orca Frank’s and Anneke Dittmer’s talks are very good.

The example here is looking a C-H bond dissociation energies, which are quite useful to study sites of metabolism. Typically I use ML to get a rough idea of which bonds to look at, then follow up with QM. The writing of multiple orca scripts was getting a bit tedious, so to facilitate the calculations I’ve used Julia to write a compound script.

julia CH_bde_calc.jl <file.xyz> <index>

There are no imports as it’s just using base module.

filename = ARGS[1]
index = parse(Int, ARGS[2])
outfilename = "radical"

# Function to extract atomic coordinates from an XYZ file
function split_xyz(filename::String)
    file_contents = readlines(filename)
    atoms = String[]
    pattern = r"([A-Za-z]{1,2})\s+([\d.-]+)\s+([\d.-]+)\s+([\d.-]+)"
    for line in file_contents
        if occursin(pattern, line)
            push!(atoms, line)
        end
    end
    return atoms
end

# Get atom list from the XYZ file
arr = split_xyz(filename)

# Remove the atom at the specified index
new_arr = copy(arr)
deleteat!(new_arr, index)

I’ve just used PBE0 as the functional for speed. Variable are defined first, then there are three jobs; calculation of energy of the whole molecule then calculations of the fragment energies. Finally we read the property files and use the three energies to calculate the BDE.

# Define functional
functional = "PBE0"

# Convert arrays to newline-separated strings
parent_xyz = join(arr, "\n")
radical_xyz = join(new_arr, "\n")

# Output filenames
basename1 = "$(outfilename)_Compound_1"
basename2 = "$(outfilename)_Compound_2"
basename3 = "$(outfilename)_Compound_3"

# Prepare script for ORCA compound mode
script = """
%Compound
  # ------------------------------------
  # First job: Get energy of parent molecule
  # ------------------------------------
  Variable HartToKcal = 627.5096080305927;  # Hartree to kcal/mol conversion factor
  Variable molec = 0.0;
  Variable frag = 0.0;
  Variable Hatom = 0.0;
  Variable basename1 = "$basename1";
  Variable basename2 = "$basename2";
  Variable basename3 = "$basename3";
  Variable energyfrag, bde, deltaG;
  Variable myProperty = "DFT_Total_en";

  New_Step
    !$functional
    ! opt anfreq
    *xyz 0 1
$parent_xyz
    *
  Step_End

  # ------------------------------------
  # Second job: Get energy of radical
  # ------------------------------------
  New_Step
    !$functional
    ! opt anfreq
    *xyz 0 2
$radical_xyz
    *
  Step_End

  # ------------------------------------
  # Third job: Get energy of hydrogen atom
  # ------------------------------------
  New_Step
    !$functional
    *xyz 0 2
    H 0.0 0.0 0.0
    *
  Step_End

  # Read calculated properties from ORCA output
  molec.ReadProperty(propertyName=myProperty, filename=basename1);
  frag.ReadProperty(propertyName=myProperty, filename=basename2);
  Hatom.ReadProperty(propertyName=myProperty, filename=basename3);

  # Compute bond dissociation energy
  energyfrag = Hatom + frag;
  deltaG = energyfrag - molec;
  bde = HartToKcal * deltaG;

  print("BDE: %.12lf kcal/mol\\n", bde);
EndRun
"""

# Write script to output file
open("$outfilename.cmp", "w") do file
    write(file, script)
end

Using tetrachloroethane as a test BDE comes out at 100.52 kcal/mol which seems about right ball park.

26th Mar 2025

compound, ORCA

QSAR with Chemprop
Chemprop is my favourite QSAR tool. It’s packed full of features, easy to use, and you can do some fairly sophisticated model building. Here are some notes on running predictions with evidential uncertainty.

Training data

I’ll use the in vogue Biogen dataset from Fang et al. 2023 which can be conveniently accessed via https://polarishub.io . As it’s a toy model I’ll assume the curated data is fine and won’t spend time finding a good data split, tailoring the descriptors, model parameters etc.
```
pip install git+https://github.com/chemprop/chemprop.git
pip install polaris-lib
pip install -U "ray[data,train,tune,serve]"
```
We’ll first generate the splits:
```
import pandas as pd
from sklearn.model_selection import train_test_split
import sys

df = pd.read_csv(sys.argv[1])

# Step 1: Split into training (80%) and remaining data (20%)
train_df, remaining_df = train_test_split(df, test_size=0.2)

# Step 2: Split remaining data into external validation (10%), and calibration (10%)
external_val_df, calibration_df = train_test_split(remaining_df, test_size=0.5)


# reset indices
train_df = train_df.reset_index(drop=True)
external_val_df = external_val_df.reset_index(drop=True)
calibration_df = calibration_df.reset_index(drop=True)

# Write the concatenated DataFrame to a CSV file
train_df.to_csv("training.csv", index=False)
external_val_df.to_csv("ext_val.csv", index=False)
calibration_df.to_csv("cal.csv",index=False)
```
Model building
```
chemprop hpopt --data-path training.csv -s "SMILES" --target-columns "AC_LOG SOLUBILITY PH 6.8 (ug/mL)" --task-type regression-evidential --search-parameter-keywords depth ffn_num_layers ffn_hidden_dim message_hidden_dim dropout --hpopt-save-dir hpoptv1 --multi-hot-atom-featurizer-mode ORGANIC -l evidential
```
Hyperparameter optimization is straightforward but time-consuming; if you’ve got the hardware, it can be sped up with --raytune-use-gpu --raytune-num-cpus 12. I found this very slow to do on an old MacBook Air M1 (>24 hours and still not complete), but Colab came to the rescue, and a free tier TPU instance worked really well (ca. 1 hr). The --multi-hot-atom-featurizer-mode flag with ORGANIC or, more recently, in v2.1.2, RIGR (Resonance Invariant Graph Representation) is more targeted and thus gives better results. At the time of writing, RIGR is a very new addition; so far, in my limited experience, the ORGANIC flag seems to slightly edge out the RIGR flag, but I expect this to be case dependent. From the pre-print, it does sound like the RIGR flag is particularly helpful for small data sets.

With the hyper-parameters identified we are ready to train the model. This can be done with the following command. Number of replicates is analogous to ‘outer loop’ of nested cross validation but at a lower cost, suitable for deep learning applications. Again if the hardware allows the training can be speeded up with --accelerator gpu. Depending on the task it important to specify the appropriate loss function (here evidential), apparently the method choice should pick this up but I am not sure it always does. With evidential regression the --evidential-regularization is model dependent, default is 0.0 and the default recommended by Soleimany et al. (2021) is 0.2
```
chemprop train --data-path training.csv -s "SMILES" --target-columns "LOG SOLUBILITY PH 6.8 (ug/mL)" --task-type regression-evidential --config-path ./hpoptv1/best_config.toml --multi-hot-atom-featurizer-mode ORGANIC --num-replicates 10 -o model_evidential --aggregation norm --epochs 50 -l evidential --evidential-regularization 0.2
```
Final step is get the predictions on the external validation set together with the uncertainty estimation which is calibrated with conformal-regression.
```
chemprop predict -s "SMILES" --test-path ext_val.csv --model-path model_evidential --preds-path preds.csv --multi-hot-atom-featurizer-mode ORGANIC --uncertainty-method evidential-total --calibration-method conformal-regression --cal-path cal.csv
```
The final results don’t look too impressive, but maybe not unsurprising; it’s a very unsophisticated model, the data range is pretty small with logS, it’s a tricky end point to predict, the experimental values can vary a lot depending on sample, the selected data isn’t well distributed over the experimental range and furthermore the applicability domain calculations suggest test set looks pretty different to train. Uncertainty estimations don’t seem to rescue the model either, or look calibrated… I’ll have to more flattering dataset next time but at least it look realistic. Next step is to see if it can be competitive TabPFN.
15th Mar 2025

chemprop, QSAR

Model Uncertainty part2: Applicability Domains

Two simple approaches I like are the sum of distance-weighted contributions (SDC) and calculating the Mahalanobis distance from the training set.

Sum of distance-weighted contributions (SDC) is the easiest to calculate;

Where d_i is the distance between the ith molecule in the training set and the test molecule, a is adjustable parameter that modulates the distance penalty (typically 3) for molecule further away from the training set. Note that d is tanimoto distance, i.e. one minus tanimoto similarity.

def calc_sdc(ref_fps, query_fps):
    d = pariwise_distances(ref_fps, query_fps, metric="jaccard",n_job=-1, force_all_finite=True)
    sdc = np.sum(np.exp(-3*d/(1-d)),axis=0)

    return sdc

Fingerprints are ideally calculated by the same method as used to build the model.

Mahalanobis distance is definitely more difficult to say, and slightly trickier to calculate. It’s simply a measure of the distance between a point and a distribution . There are good implications for in python, Julia , etc. but something I thought would be useful to do would be to use the Facebook A.I. similarity search (FAISS) tool as it is super fast, easy to use, and so enabling application to large libraries.

def create_vector(vector):
    """
    Transform a vector to follow a unit Gaussian distribution using Mahalanobis transformation.
    """
    vector = np.ascontiguousarray(vector)
    xc = vector - vector.mean(0)
    cov = np.dot(xc.T, xc) / xc.shape[0]
    cov += 2 * np.eye(cov.shape[0])
    L = np.linalg.cholesky(cov)
    mahalanobis_transform = np.linalg.inv(L)
    return np.dot(vector, mahalanobis_transform.T)

With the vectors created from this function we create an in memory database of the transformed training vectors which we will then search.

def create_index(input_vector):
    """
    Create a FAISS index for a given set of vectors.
    """
    vector_dimension = input_vector.shape[1]
    index = faiss.IndexFlatL2(vector_dimension)
    index.add(input_vector)
    print(f"Index created with {index.ntotal} vectors")
    return index

Finally we search the indexed database and write the distances to file.

def main():
    if len(sys.argv) != 3:
        print("Usage: python AD_calc_mahalanobis.py <training_csv> <predict_csv>")
        sys.exit(1)

    training_df = pd.read_csv(sys.argv[1])
    predict_df = pd.read_csv(sys.argv[2])

    training_smiles = training_df["smiles"].tolist()
    predict_smiles = predict_df["smiles"].tolist()

    train_fps = smiles_to_fingerprint_array(training_smiles)
    predict_fps = smiles_to_fingerprint_array(predict_smiles)
    ###### SDC #############################################################


    sdc = calc_sdc(predict_fps, train_fps)

    predict_df["scd"] = sdc


    ######## mahalanobisDistance #################################

    array_train = create_vector(train_fps)
    array_predict = create_vector(predict_fps)


    index = create_index(array_train)
    distances, _ = index.search(array_predict, 1)

    predict_df["mahalanobisDistance"] = distances.flatten()

    ############# Write output ####################################
    predict_df.to_csv("test_cmps_with_distance.csv", index=False)
    print("Results saved to test_cmps_with_distance.csv")

if __name__ == "__main__":
    main()

To sanity check everything we will plot the SDC of the external test set, calibration set and a 10000 random smiles generated using MoLeR sampling a chembl derived model provided. The box-plot was generate with AlgebraOfGraphics.

using Colors
using AlgebraOfGraphics, CairoMakie
using DataFrames, CSV

update_theme!(
           fontsize=16,
           markersize=8,
           Axis=(title=L"Box Plot SDC",))

# define colour scheme 
Sun_Orange_Bright = RGB(252/255, 229/255, 0)
Plant_Green_Bright = RGB(80/255, 230/255, 0)
Air_Blue_Bright = RGB(0, 230/255, 190/255)
Energy_Pink_Bright = RGB(255/255, 0, 100/255)

#read data
df_1 = CSV.read(ARGS[1], DataFrame)
df_2 = CSV.read(ARGS[2], DataFrame)
df_3 = CSV.read(ARGS[3], DataFrame)


y1 = df_1[!, :scd]
y2 = df_2[!, :scd]
y3 = df_3[!, :scd]

x = vcat(fill("test", length(y1)), fill("cal", length(y2)), fill("random", length(y3)))  # Group labels
y = vcat(y1, y2, y3)  # Combine columns into one vector

# Create a named tuple for AlgebraOfGraphics
df = (; x, y)
plt = data(df) * mapping(:x=>L"x", :y =>L"y", color=:x) * visual(BoxPlot) * data(df) * mapping(:x=>L"x", :y =>L"y", color=:x) * visual(BoxPlot)

colors = [Plant_Green_Bright, Air_Blue_Bright, Energy_Pink_Bright]

fg = draw(plt, scales(Color = (; palette = colors)))

save("BoxPlot.svg",fg, px_per_unit = 3)

As expected test set and calibration set look similar, whereas the random smiles looks more different.

15th Feb 2025

AD, AlgebraOfGraphics, BoxPlot, SDC, uncertainty

Extracting uncertainty estimations from raw data.

I recently came across this blog article from Dr. Lewis Martin on Bayesian IC₅₀ fitting. The determination of uncertainty from raw data is something that is of interest to me. In the article Lewis use PyMC3, which is pretty popular, the analogues package in Julia is Turing so we will use that.

Imports

StaticArray.jl is mainly there to speed up the code, there maybe be better ways of doing this with threading or using a GPU.

using AlgebraOfGraphics, CairoMakie
using DataFrames
using  Turing, Distributions
using StaticArrays

Generate Fake Data

First define the concentration range for the x-axis.

# x points are in units of µM - they describe the 
# test ligand concentrations across a 9pt CRC
concs = [0.01, 0.032, 0.1, 0.32, 1, 3.2, 10, 32, 100]
# Convert concentrations to log10 Molar units
log_concs = -log10.(concs ./ 1e6)

Make up some pIC50 data for the y-axis using a sigmoid function.

function sigmoid(x; middle, slope, top, bottom)
    return top / (1 + ℯ^(-(x - middle) * slope) - bottom)
end

sigmoid_values = sigmoid.(log_concs; middle=6.1, slope=2, top=100, bottom=0)

We’ll check the data by plotting the data, here we will use Algebra of Graphics

# plot initial data.
df = (; log_concs, sigmoid_values, )
xy = data(df) * mapping(:log_concs => L"concentrations", :sigmoid_values => L"Activity")
layers = visual(Scatter)
fg = draw(layers * xy)
save("figure1.png", fg, px_per_unit = 3)

Next we’ll re-plot the data with duplicates and noise.

#add some duplicates and noise.
concs = repeat(log_concs, inner=2)        
responses = repeat(sigmoid_values, inner=2)  
# Add random noise (Normal distribution with mean 0, std dev 5)
responses .+= randn(length(responses)) .* 5

With the data in hand we ready for modelling, first we define a model. Priors for top and bottom are inherently set to 100 and 0. For the variance we use the inverse gamma, the alternative would be the truncated normal err ~ Truncated(Normal(0, 20), 0, Inf). After this we build the model and then do a “No U–Turn Sampling” with 10,000 datapoints.

# Define the Turing model
@model function sigmoid_model(concs, responses)
    concs = SVector{length(concs)}(concs)
    responses = SVector{length(responses)}(responses)
    # Priors
    s² ~ InverseGamma(0.5, 20)
    midpoint ~ Normal(5, sqrt(s²))
    slope ~ Normal(0, sqrt(s²))
    bottom ~ Normal(0, sqrt(s²))
    top ~ Normal(100, sqrt(s²))
    
    # Sigmoid function: 
    proportion = @. top / (1 + exp(-(concs + midpoint) * slope)) - bottom
    # Likelihood (observed data with Gaussian noise)
    responses ~ MvNormal(proportion, s²)
end
end

model = sigmoid_model(concs, responses)
chain = sample(model, NUTS(), 10_000; tune=3000, target_accept=0.9)

Finally plot the data as density plot, which looks good mean is approximately 6.1 and standard deviation is approximately 0.1.

df = DataFrame(chain)
df.midpoint = abs.(df.midpoint)
xy = data(df) * mapping(:midpoint => L"midpoint")
layers=AlgebraOfGraphics.density()
fg = draw(layers * xy)
save("figure3.png", fg, px_per_unit = 3)

1st Jan 2025

Bayesian, Julia, Turing