<jupyter_start><jupyter_text>Scaffold splitting and initial model training using chemprop> Splitting data using a chemical space aware technique - scaffold splitting. Data is then trained using a directed message passing neural network (D-MPNN).- toc: true- badges: true- comments: true- categories: [python, bioinformatics, datasets, SMILES, cheminformatics, datamol, RDKit, molecules]- image: images/play.gif  ChempropI found the chemprop library in the paper [A Deep Learning Approach to Antibiotic Discovery]("https://www.cell.com/cell/fulltext/S0092-8674(20)30102-1"), it uses the directed message passing neural network (D-MPNN) described in [Analyzing Learned Molecular Representations for Property Prediction]("https://pubs.acs.org/doi/abs/10.1021/acs.jcim.9b00237").The drug discovery workflow described in the paper uses a training dataset of 2335 molecules with binary e.coli growth inhibition labels to train a classification model.The D-MPNN architecture as described in the paper:> The D-MPNN architecture translates the graph representation of a molecule into a continuous vector via a directed bond-based message-passing approach. This builds a molecular representation by iteratively aggregating the features of individual atoms and bonds. The model operates by passing “messages” along bonds that encode information about neighboring atoms and bonds. By applying this message passing operation multiple times, the model constructs higher-level bond messages that contain information about larger chemical substructures. The highest-level bond messages are then combined into a single continuous vector representing the entire molecule.The researchers generated molecular features using RDKit to tackle overfitting, increased algorithm robustness by using an ensemble of classifiers, and optimized hyperparameters with bayesian optimization. The model achieved a ROC-AUC of 0.896 on test data.This notebook will apply the neural network to the dataset of blood-brain barrier molecules gathered in the previous notebooks.  Scaffold splittingA distributional shift is a change in the data distribution between an algorithm's training dataset, and a dataset it encounters when deployed. Distributional shifts are common in the field of molecular property prediction where chemical space is huge and different areas of it exhibit high structural heterogeneity. The shifts tend to be large and difficult to handle for machine learning models. Scaffold splitting is one solution to this distributional shift. It's first described in ["The properties of known drugs. 1. Molecular frameworks]("https://pubs.acs.org/doi/10.1021/jm9602928"):  > A molecular scaffold reduces the chemical structure of a compound to its core components, essentially by removing all side chains and only keeping ring systems and parts that link together ring systems. An additional option for making molecular scaffolds even more general is to “forget” the identities of the bonds and atoms by replacing all atoms with carbons and all bonds with single bonds.A Murcko scaffold essentially extracts the molecular backbone by extracting the ring structures and the linkers that connect them. These scaffolds collapse molecules into bins based on their generated Murcko scaffold, reducing the number of highly similar structures. When splitting the data by scaffold, molecules sharing a scaffold are in the same split, meaning no similar molecules will be found across different splits. Any bins larger than half of the desired test set size are placed into the training set, to guarantee the scaffold diversity of the validation and test sets. All remaining bins are placed randomly into the training, validation, and test sets until each set has reached its desired size. Clustering molecules by scaffold ensures that the model will be trained on structurally diverse folds when cross-validating.Compared to a random split, a scaffold split is a more challenging and realistic evaluation setting as it more closely approximates the split present in real-world property prediction data.A molecule compared to its scaffold below:  Molecular compound and its generic Bemis-Murcko scaffold<jupyter_code>from rdkit import Chem
from rdkit.Chem.Scaffolds import MurckoScaffold

# define compound via its SMILES string
smiles = "CN1CCCCC1CCN2C3=CC=CC=C3SC4=C2C=C(C=C4)SC"
# convert SMILES string to RDKit mol object 
mol = Chem.MolFromSmiles(smiles)
# create RDKit mol object corresponding to Bemis-Murcko scaffold of original compound
mol_scaffold = MurckoScaffold.GetScaffoldForMol(mol)
# make the scaffold generic by replacing all atoms with carbons and all bonds with single bonds
mol_scaffold_generic = MurckoScaffold.MakeScaffoldGeneric(mol_scaffold)
# convert the generic scaffold mol object back to a SMILES string format
smiles_scaffold_generic = Chem.CanonSmiles(Chem.MolToSmiles(mol_scaffold_generic))
# display compound and its generic Bemis-Murcko scaffold
display(mol)
print(smiles)
display(mol_scaffold_generic)
print(smiles_scaffold_generic)<jupyter_output><empty_output><jupyter_text>Generating generic Bemis-Murcko scaffoldsBelow I've written code to generate the scaffolds and a few other related molecular representations. But binning and writing the code to split them myself would take too much time so I use a built-in feature of chemprop to split the data before model training.<jupyter_code># import packages

import chemprop
import pandas as pd
import datamol as dm
pd.options.mode.chained_assignment = None  # default='warn'


# loading the datasets
MolNet = pd.read_csv("data/MoleculeNet.csv")
B3DB = pd.read_csv("data/B3DB.csv")
b3_molecules = pd.read_csv("data/b3_molecules.csv")<jupyter_output><empty_output><jupyter_text>scaffold_split functionThis function generates and sanitizes mols, generates a scaffold from the mol object, generalizes the scaffold, and lastly converts the scaffold to a SMILES string, returning a dataframe.<jupyter_code>def scaffold_split(df):
    df["mol"] = [Chem.MolFromSmiles(x) for x in df["standard_smiles"]] # generating mols from the standard_smiles column
    df["mol"] = [dm.sanitize_mol(x, sanifix=True, charge_neutral=True) for x in df['mol']] # sanitize mol objects
    df = df.dropna() # dropping NA values
    df["scaffold"] = [MurckoScaffold.GetScaffoldForMol(x) for x in df["mol"]] # generating scaffolds from mol object
    df["mol_scaffold_generic"] = [MurckoScaffold.MakeScaffoldGeneric(x) for x in df["scaffold"]] # generalizing scaffolds
    # convert the generic scaffold mol object back to a SMILES string format
    df["smiles_scaffold_generic"] = [Chem.CanonSmiles(Chem.MolToSmiles(x)) for x in df["mol_scaffold_generic"]]
    
    return df

# the data prior to processing
MolNet.head(1)

# results of the scaffold generating function
data_split = scaffold_split(b3_molecules)
data_split.head(1)<jupyter_output><empty_output><jupyter_text>chemprop_prep functionThis function serves to prepare the data to be used by the chemprop library. Specifically, the input data must contain a SMILES string in the first column and the target value - a binary value in this experiment. The resulting dataframe is saved as a CSV file to later be used by chemprop.<jupyter_code># A function to process the data for use with the chemprop library
def chemprop_prep(df, filename):
    df = df.drop(["mol", "SMILES", "selfies", "inchi", "inchikey"], axis=1) # drop all columns except the smiles and target
    df["smiles"] = df["standard_smiles"] # use standard smiles inplace of smiles
    df = df.drop(["standard_smiles"], axis=1) # drop this column now
    df = df[["smiles", "BBB+/BBB-"]] # reorder the columns with smiles first and target second
    df.to_csv('./data/' + filename + '.csv', index=False) # save the file

    return df

# Processing the three different datasets
molnet_chemprop = chemprop_prep(MolNet, 'molnet_chemprop')
B3DB_chemprop = chemprop_prep(B3DB, 'B3DB_chemprep')
b3_mol_chemprop = chemprop_prep(b3_molecules, 'b3_mol_chemprop')

# results of processing
molnet_chemprop.head(1)<jupyter_output><empty_output><jupyter_text>Model training with chempropI ran three experiments with slight differences to test out the predictive ability on the blood-brain barrier molecules:1. The split type is set to scaffold balanced; the number of folds is 10, class balancing is set to true, and 200 features are generated with RDKits 2d descriptor generator. The features are scaled when calculated and thus no further scaling is specified.2. This run sets the number of training epochs within each fold at 100. The rest of the settings mentioned above are kept the same.3. The last run has the number of folds increased from 10 to 30 while epochs per fold are set to 20.These models are trained using the MoleculeNet data (~2035 molecules), further tests will use the combined dataset with ~8100 molecules.  Model training 1<jupyter_code>arguments = [
    '--data_path', './data/molnet_chemprop.csv',
    '--dataset_type', 'classification',
    '--save_dir', './data/chemprop_checkpoints/',
    '--split_type', 'scaffold_balanced',
 #   '--separate_val_path', './data/chemprop_B3DB.csv',
    '--num_folds', '10',
    '--class_balance',
    '--features_generator', 'rdkit_2d_normalized',
    '--no_features_scaling',
    '--quiet'
    

]

args = chemprop.args.TrainArgs().parse_args(arguments)
mean_score, std_score = chemprop.train.cross_validate(args=args, train_func=chemprop.train.run_training)<jupyter_output><empty_output><jupyter_text>Model training 2<jupyter_code>arguments = [
    '--data_path', './data/molnet_chemprop.csv',
    '--dataset_type', 'classification',
    '--save_dir', './data/chemprop_checkpoints/',
    '--split_type', 'scaffold_balanced',
 #   '--separate_val_path', './data/chemprop_B3DB.csv',
    '--num_folds', '10',
    '--class_balance',
    '--features_generator', 'rdkit_2d_normalized',
    '--no_features_scaling',
    '--quiet',
    '--epochs', '100'
    

]

args = chemprop.args.TrainArgs().parse_args(arguments)
mean_score, std_score = chemprop.train.cross_validate(args=args, train_func=chemprop.train.run_training)<jupyter_output><empty_output><jupyter_text>Model training 3<jupyter_code>arguments = [
    '--data_path', './data/molnet_chemprop.csv',
    '--dataset_type', 'classification',
    '--save_dir', './data/chemprop_checkpoints/',
    '--split_type', 'scaffold_balanced',
 #   '--separate_val_path', './data/chemprop_B3DB.csv',
    '--num_folds', '30',
    '--class_balance',
    '--features_generator', 'rdkit_2d_normalized',
    '--no_features_scaling',
    '--quiet',
    '--epochs', '20',
    

]

args = chemprop.args.TrainArgs().parse_args(arguments)
mean_score, std_score = chemprop.train.cross_validate(args=args, train_func=chemprop.train.run_training)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07905943343825217, 'fraction_numerical': 0.011109053878911312, 'mean_word_length': 5.269878119558909, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 32, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '138831', 'license': 'Apache-2.0', 'n_tokens_mistral': 3237, 'n_tokens_neox': 2943, 'n_words': 1173, 'path': '_notebooks/2021-11-12-data-splitting-and-training.ipynb.ipynb', 'repo_name': 'gdeol4/neongenes', 'size': '10802'}