<jupyter_start><jupyter_text>T007 · Ligand-based screening: machine learningAuthors:* <NAME>, CADD seminar 2018, Charité/FU Berlin* <NAME>, CADD seminar 2018, Charité/FU Berlin* <NAME>, 2019-2020, [Volkamer lab](https://volkamerlab.org)* <NAME>, 2019-2020, [Volkamer lab](https://volkamerlab.org) __Talktorial T007__: This talktorial is part of the TeachOpenCADD pipeline described in the [first TeachOpenCADD paper](https://jcheminf.biomedcentral.com/articles/10.1186/s13321-019-0351-x), comprising of talktorials T001-T010.  Aim of this talktorialDue to larger available data sources, machine learning (ML) gained momentum in drug discovery and especially in ligand-based virtual screening. In this talktorial, we learn how to use different supervised ML algorithms to predict the activity of novel compounds against our target of interest (EGFR).  Contents in _Theory_* Data preparation: Molecule encoding* Machine learning (ML)    * Supervised learning* Model validation and evaluation    * Validation strategy: K-fold cross-validation    *  Performance measures  Contents in _Practical_* Load compound and activity data* Data preparation    * Data labeling    * Molecule encoding* Machine learning    * Helper functions    * Random forest classifier    * Support vector classifier    * Neural network classifier    * Cross-validation  References* "Fingerprints in the RDKit" [slides](https://www.rdkit.org/UGM/2012/Landrum_RDKit_UGM.Fingerprints.Final.pptx.pdf), <NAME>, RDKit UGM 2012* Extended-connectivity fingerprints (ECFPs): Rogers, David, and <NAME>. "Extended-connectivity fingerprints." [_Journal of chemical information and modeling_ 50.5 (2010): 742-754.](https://doi.org/10.1021/ci100050t)* Machine learning (ML):  * Random forest (RF): <NAME>. "Random Forests". [_Machine Learning_ **45**, 5–32 (2001).](https://link.springer.com/article/10.1023%2FA%3A1010933404324)  * Support vector machines (SVM): <NAME>., <NAME>. "Support-vector networks". [_Machine Learning_ **20**, 273–297 (1995).](https://link.springer.com/article/10.1007%2FBF00994018)  * Artificial neural networks (ANN): <NAME>, Marcel, and <NAME>. "Artificial neural networks as models of neural information processing." [_Frontiers in Computational Neuroscience_ 11 (2017): 114.](https://doi.org/10.3389/fncom.2017.00114)* Performance:   * Sensitivity and specificity ([Wikipedia](https://en.wikipedia.org/wiki/Sensitivity_and_specificity))  * ROC curve and AUC ([Wikipedia](https://en.wikipedia.org/wiki/Receiver_operating_characteristicArea_under_the_curve))* See also [github notebook by <NAME>](https://github.com/Team-SKI/Publications/tree/master/Profiling_prediction_of_kinase_inhibitors) from [*J. Med. Chem.*, 2017, 60, 474−485](https://pubs.acs.org/doi/10.1021/acs.jmedchem.6b01611) * Activity cutoff $pIC_{50} = 6.3$ used in this talktorial  * Profiling Prediction of Kinase Inhibitors: Toward the Virtual Assay [J. Med. Chem. (2017), 60, 474-485](https://doi.org/10.1021/acs.jmedchem.6b01611)  * Notebook accompanying the publication mentioned before: [Notebook](https://github.com/Team-SKI/Publications/blob/master/Profiling_prediction_of_kinase_inhibitors/Build_ABL1_model.ipynb)  TheoryTo successfully apply ML, we need a large data set of molecules, a molecular encoding, a label per molecule in the data set, and a ML algorithm to train a model. Then, we can make predictions for new molecules.![ML overview](images/ML_overview.png)_Figure 1_: Machine learning overview: Molecular encoding, label, ML algorithm, prediction. Figure by <NAME>.  Data preparation: Molecule encodingFor ML, molecules need to be converted into a list of features. Often molecular fingerprints are used as representation. The fingerprints used in this talktorial as implemented in RDKit (more info can be found in a  [presentation by <NAME>](https://www.rdkit.org/UGM/2012/Landrum_RDKit_UGM.Fingerprints.Final.pptx.pdf)) are:* **maccs**: 'MACCS keys are 166 bit structural key descriptors in which each bit is associated with a SMARTS pattern.' (see OpenEye's `MACCS` [docs](https://docs.eyesopen.com/toolkits/python/graphsimtk/fingerprint.htmlmaccs))* **Morgan fingerprints** (and **ECFP**): 'Extended-Connectivity Fingerprints (ECFPs) are circular topological fingerprints designed for molecular characterization, similarity searching, and structure-activity modeling.' (see ChemAxon's `ECFP` [docs](https://docs.chemaxon.com/display/docs/Extended+Connectivity+Fingerprint+ECFP)) The original implementation of the ECFPs was done in Pipeline Pilot which is not open-source. Instead we use the implementation from RDKit which is called Morgan fingerprint. The two most important parameters of these fingerprints are the radius and fingerprint length. The first specifies the radius of circular neighborhoods considered for each atom. Here two radii are considered: 2 and 3. The length parameter specifies the length to which the bit string representation is hashed. The default length is 2048.  Machine learning (ML)ML can be applied for (text adapted from [scikit-learn page](http://scikit-learn.org/stable/)):* **Classification (supervised)**: Identify which category an object belongs to (e.g. : Nearest neighbors, Naive Bayes, RF, SVM, ...)* Regression: Prediction of a continuous-values attribute associated with an object* Clustering (unsupervised): Automated grouping of similar objects into sets (see also **Talktorial T005**)  Supervised learningA learning algorithm creates rules by finding patterns in the training data. * **Random Forest (RF)**: Ensemble of decision trees. A single decision tree splits the features of the input vector in a way that maximizes an objective function. In the random forest algorithm, the trees that are grown are de-correlated because the choice of features for the splits are chosen randomly.* **Support Vector Machines (SVMs)**: SVMs can efficiently perform a non-linear classification using what is called the kernel trick, implicitly mapping their inputs into high-dimensional feature spaces. The classifier is based on the idea of maximizing the margin as the objective function.  * **Artificial neural networks (ANNs)**: An ANN is based on a collection of connected units or nodes called artificial neurons which loosely model the neurons in a biological brain. Each connection, like the synapses in a biological brain, can transmit a signal from one artificial neuron to another. An artificial neuron that receives a signal can process it and then signal additional artificial neurons connected to it.![ANN_example](images/ANN_wiki.png)_Figure 2_: Example of a neural network with one hidden layer. Figure taken from [Wikipedia](https://en.wikipedia.org/wiki/Artificial_neural_network).  Model validation and evaluation  Validation strategy: K-fold cross validation* This model validation technique splits the dataset in two groups in an iterative manner:    * Training data set: Considered as the known dataset on which the model is trained    * Test dataset: Unknown dataset on which the model is then tested    * Process is repeated k-times* The goal is to test the ability of the model to predict data which it has never seen before in order to flag problems known as over-fitting and to assess the generalization ability of the model.  Performance measures* **Sensitivity**, also true positive rate    * TPR = TP/(FN + TP)     * _Intuitively_: Out of all actual positives, how many were predicted as positive?* **Specificity**, also true negative rate    * TNR = TN/(FP + TN)    * _Intuitively_: Out of all actual negatives, how many were predicted as negative?* **Accuracy**, also the trueness    * ACC = (TP + TN)/(TP + TN + FP + FN)    * _Intuitively_: Proportion of correct predictions.* **ROC-curve**, receiver operating characteristic curve    * A graphical plot that illustrates the diagnostic ability of our classifier    * Plots the sensitivity against the specificity* **AUC**, the area under the ROC curve (AUC):      * Describes the probability that a classifier will rank a randomly chosen positive instance higher than a negative one    * Values between 0 and 1, the higher the better | What the model predicts  | True active  |  True inactive ||---|---|---|| active  |  True Positive (TP) |  False Positive (FP) || inactive  |  False Negative (FN) |  True Negative (TN) |  Practical<jupyter_code>from pathlib import Path
from warnings import filterwarnings
import time

import pandas as pd
import numpy as np
from sklearn import svm, metrics, clone
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import RandomForestRegressor
from sklearn.neural_network import MLPClassifier
from sklearn.model_selection import KFold, train_test_split
from sklearn.metrics import auc, accuracy_score, recall_score
from sklearn.metrics import roc_curve, roc_auc_score
import matplotlib.pyplot as plt
from rdkit import Chem
from rdkit.Chem import MACCSkeys
from rdkit.Chem.AllChem import GetMorganFingerprintAsBitVect

from teachopencadd.utils import seed_everything

# Silence some expected warnings
filterwarnings("ignore")
# Fix seed for reproducible results
SEED = 22
seed_everything(SEED)
# Set path to this notebook
HERE = Path(_dh[-1])
DATA = HERE / "data"<jupyter_output><empty_output><jupyter_text>Load compound and activity dataLet's start by loading our data, which focuses on the Epidermal growth factor receptor (EGFR) kinase. The *csv* file from **Talktorial T002** is loaded into a dataframe with the important columns:* CHEMBL-ID* SMILES string of the corresponding compound* Measured affinity: pIC50<jupyter_code># Read data from previous talktorials
chembl_df = pd.read_csv(
    HERE / "../T002_compound_adme/data/A2A_compounds_lipinski.csv",
    index_col=0,
)

# Look at head
print("Shape of dataframe : ", chembl_df.shape)
chembl_df.head()
# NBVAL_CHECK_OUTPUT
# Keep only the columns we want
chembl_df = chembl_df[["molecule_chembl_id", "smiles", "pKi"]]
chembl_df.head()
# NBVAL_CHECK_OUTPUT<jupyter_output><empty_output><jupyter_text>Data preparation  Data labelingWe need to classify each compound as active or inactive. Therefore, we use the pIC50 value.* pIC50 = -log10(IC50) * IC50 describes the amount of substance needed to inhibit, _in vitro_, a process by 50% .* A common cut-off value to discretize pIC50 data is 6.3, which we will use for our experiment (refer to [<NAME>. Chem. (2017), 60, 474-485](https://doi.org/10.1021/acs.jmedchem.6b01611) and the corresponding[notebook](https://github.com/Team-SKI/Publications/blob/master/Profiling_prediction_of_kinase_inhibitors/Build_ABL1_model.ipynb))* Note that there are several other suggestions for an activity cut-off ranging from an pIC50 value of 5 to 7 in the literature or even to define an exclusion range when not to take data points.<jupyter_code># Add column for activity
chembl_df["active"] = np.zeros(len(chembl_df))

# Mark every molecule as active with an pKi of >= 6.3, 0 otherwise
chembl_df.loc[chembl_df[chembl_df.pKi >= 6.3].index, "active"] = 1.0

# NBVAL_CHECK_OUTPUT
print("Number of active compounds:", int(chembl_df.active.sum()))
print("Number of inactive compounds:", len(chembl_df) - int(chembl_df.active.sum()))
chembl_df.head()
# NBVAL_CHECK_OUTPUT<jupyter_output><empty_output><jupyter_text>Molecule encodingNow we define a function `smiles_to_fp` to generate fingerprints from SMILES.For now, we incorporated the choice between the following fingerprints:* maccs* morgan2 and morgan3<jupyter_code>def smiles_to_fp(smiles, method="maccs", n_bits=2048):
    """
    Encode a molecule from a SMILES string into a fingerprint.

    Parameters
    ----------
    smiles : str
        The SMILES string defining the molecule.

    method : str
        The type of fingerprint to use. Default is MACCS keys.

    n_bits : int
        The length of the fingerprint.

    Returns
    -------
    array
        The fingerprint array.

    """

    # convert smiles to RDKit mol object
    mol = Chem.MolFromSmiles(smiles)

    if method == "maccs":
        return np.array(MACCSkeys.GenMACCSKeys(mol))
    if method == "morgan2":
        return np.array(GetMorganFingerprintAsBitVect(mol, 2, nBits=n_bits))
    if method == "morgan3":
        return np.array(GetMorganFingerprintAsBitVect(mol, 3, nBits=n_bits))
    else:
        # NBVAL_CHECK_OUTPUT
        print(f"Warning: Wrong method specified: {method}. Default will be used instead.")
        return np.array(MACCSkeys.GenMACCSKeys(mol))
compound_df = chembl_df.copy()
# Add column for fingerprint
compound_df["fp"] = compound_df["smiles"].apply(smiles_to_fp)
compound_df.head(3)
# NBVAL_CHECK_OUTPUT<jupyter_output><empty_output><jupyter_text>Machine Learning (ML)In the following, we will try several ML approaches to classify our molecules. We will use:* Random Forest (RF)* Support Vector Machine (SVM) * Artificial Neural Network (ANN) Additionally, we will comment on the results.The goal is to test the ability of the model to predict data which it has never seen before in order to flag problems known as over fitting and to assess the generalization ability of the model.We start by defining a function `model_training_and_validation` which fits a model on a random train-test split of the data and returns measures such as accuracy, sensitivity, specificity and AUC evaluated on the test set. We also plot the ROC curves using `plot_roc_curves_for_models`.We then define a function named `crossvalidation` which executes a cross validation procedure and prints the statistics of the results over the folds.  Helper functionsHelper function to plot customized ROC curves. Code inspired by [stackoverflow](https://stackoverflow.com/questions/42894871/how-to-plot-multiple-roc-curves-in-one-plot-with-legend-and-auc-scores-in-python).<jupyter_code>def plot_roc_curves_for_models(models, test_x, test_y, save_png=False):
    """
    Helper function to plot customized roc curve.

    Parameters
    ----------
    models: dict
        Dictionary of pretrained machine learning models.
    test_x: list
        Molecular fingerprints for test set.
    test_y: list
        Associated activity labels for test set.
    save_png: bool
        Save image to disk (default = False)

    Returns
    -------
    fig:
        Figure.
    """

    fig, ax = plt.subplots()

    # Below for loop iterates through your models list
    for model in models:
        # Select the model
        ml_model = model["model"]
        # Prediction probability on test set
        test_prob = ml_model.predict_proba(test_x)[:, 1]
        # Prediction class on test set
        test_pred = ml_model.predict(test_x)
        # Compute False postive rate and True positive rate
        fpr, tpr, thresholds = metrics.roc_curve(test_y, test_prob)
        # Calculate Area under the curve to display on the plot
        auc = roc_auc_score(test_y, test_prob)
        # Plot the computed values
        ax.plot(fpr, tpr, label=(f"{model['label']} AUC area = {auc:.2f}"))

    # Custom settings for the plot
    ax.plot([0, 1], [0, 1], "r--")
    ax.set_xlabel("False Positive Rate")
    ax.set_ylabel("True Positive Rate")
    ax.set_title("Receiver Operating Characteristic")
    ax.legend(loc="lower right")
    # Save plot
    if save_png:
        fig.savefig(f"{DATA}/roc_auc", dpi=300, bbox_inches="tight", transparent=True)
    return fig<jupyter_output><empty_output><jupyter_text>Helper function to calculate model performance.<jupyter_code>def model_performance(ml_model, test_x, test_y, verbose=True):
    """
    Helper function to calculate model performance

    Parameters
    ----------
    ml_model: sklearn model object
        The machine learning model to train.
    test_x: list
        Molecular fingerprints for test set.
    test_y: list
        Associated activity labels for test set.
    verbose: bool
        Print performance measure (default = True)

    Returns
    -------
    tuple:
        Accuracy, sensitivity, specificity, auc on test set.
    """

    # Prediction probability on test set
    test_prob = ml_model.predict_proba(test_x)[:, 1]

    # Prediction class on test set
    test_pred = ml_model.predict(test_x)

    # Performance of model on test set
    accuracy = accuracy_score(test_y, test_pred)
    sens = recall_score(test_y, test_pred)
    spec = recall_score(test_y, test_pred, pos_label=0)
    auc = roc_auc_score(test_y, test_prob)

    if verbose:
        # Print performance results
        # NBVAL_CHECK_OUTPUT        print(f"Accuracy: {accuracy:.2}")
        print(f"Sensitivity: {sens:.2f}")
        print(f"Specificity: {spec:.2f}")
        print(f"AUC: {auc:.2f}")

    return accuracy, sens, spec, auc<jupyter_output><empty_output><jupyter_text>Helper function to fit a machine learning model on a random train-test split of the data and return the performance measures.<jupyter_code>def model_training_and_validation(ml_model, name, splits, verbose=True):
    """
    Fit a machine learning model on a random train-test split of the data
    and return the performance measures.

    Parameters
    ----------
    ml_model: sklearn model object
        The machine learning model to train.
    name: str
        Name of machine learning algorithm: RF, SVM, ANN
    splits: list
        List of desciptor and label data: train_x, test_x, train_y, test_y.
    verbose: bool
        Print performance info (default = True)

    Returns
    -------
    tuple:
        Accuracy, sensitivity, specificity, auc on test set.

    """
    train_x, test_x, train_y, test_y = splits

    # Fit the model
    ml_model.fit(train_x, train_y)

    # Calculate model performance results
    accuracy, sens, spec, auc = model_performance(ml_model, test_x, test_y, verbose)

    return accuracy, sens, spec, auc<jupyter_output><empty_output><jupyter_text>**Preprocessing**: Split the data (will be reused for the other models)<jupyter_code>fingerprint_to_model = compound_df.fp.tolist()
label_to_model = compound_df.active.tolist()

# Split data randomly in train and test set
# note that we use test/train_x for the respective fingerprint splits
# and test/train_y for the respective label splits
(
    static_train_x,
    static_test_x,
    static_train_y,
    static_test_y,
) = train_test_split(fingerprint_to_model, label_to_model, test_size=0.2, random_state=SEED)
splits = [static_train_x, static_test_x, static_train_y, static_test_y]
# NBVAL_CHECK_OUTPUT
print("Training data size:", len(static_train_x))
print("Test data size:", len(static_test_x))<jupyter_output>Training data size: 2836
Test data size: 709<jupyter_text>Random forest classifierWe start with a random forest classifier, where we first set the parameters. We train the model on a random train-test split and plot the results.<jupyter_code># Set model parameter for random forest
param = {
    "n_estimators": 100,  # number of trees to grows
    "criterion": "entropy",  # cost function to be optimized for a split
}
model_RF = RandomForestClassifier(**param)
# Fit model on single split
performance_measures = model_training_and_validation(model_RF, "RF", splits)
# Initialize the list that stores all models. First one is RF.
models = [{"label": "Model_RF", "model": model_RF}]
# Plot roc curve
plot_roc_curves_for_models(models, static_test_x, static_test_y);<jupyter_output><empty_output><jupyter_text>Support vector classifierHere we train a SVM with a radial-basis function kernel (also: squared-exponential kernel). For more information, see [sklearn RBF kernel](http://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.kernels.RBF.html).<jupyter_code># Specify model
model_SVM = svm.SVC(kernel="rbf", C=1, gamma=0.1, probability=True)

# Fit model on single split
performance_measures = model_training_and_validation(model_SVM, "SVM", splits)
# Append SVM model
models.append({"label": "Model_SVM", "model": model_SVM})
# Plot roc curve
plot_roc_curves_for_models(models, static_test_x, static_test_y);<jupyter_output><empty_output><jupyter_text>Neural network classifierThe last approach we try here is a neural network model. We train an MLPClassifier (Multi-layer Perceptron classifier) with 3 layers, each with 5 neurons. As before, we do the crossvalidation procedure and plot the results. For more information on MLP, see [sklearn MLPClassifier](http://scikit-learn.org/stable/modules/generated/sklearn.neural_network.MLPClassifier.html).<jupyter_code># Specify model
model_ANN = MLPClassifier(hidden_layer_sizes=(5, 3), random_state=SEED)

# Fit model on single split
performance_measures = model_training_and_validation(model_ANN, "ANN", splits)
# Append ANN model
models.append({"label": "Model_ANN", "model": model_ANN})
# Plot roc curve
plot_roc_curves_for_models(models, static_test_x, static_test_y, True);<jupyter_output><empty_output><jupyter_text>Our models show very good values for all measured values (see AUCs) and thus seem to be predictive.  Cross-validationNext, we will perform cross-validation experiments with the three different models.Therefore, we define a helper function for machine learning model training and validation in a cross-validation loop.<jupyter_code>def crossvalidation(ml_model, df, n_folds=5, verbose=False):
    """
    Machine learning model training and validation in a cross-validation loop.

    Parameters
    ----------
    ml_model: sklearn model object
        The machine learning model to train.
    df: pd.DataFrame
        Data set with SMILES and their associated activity labels.
    n_folds: int, optional
        Number of folds for cross-validation.
    verbose: bool, optional
        Performance measures are printed.

    Returns
    -------
    None

    """
    t0 = time.time()
    # Shuffle the indices for the k-fold cross-validation
    kf = KFold(n_splits=n_folds, shuffle=True, random_state=SEED)

    # Results for each of the cross-validation folds
    acc_per_fold = []
    sens_per_fold = []
    spec_per_fold = []
    auc_per_fold = []

    # Loop over the folds
    for train_index, test_index in kf.split(df):
        # clone model -- we want a fresh copy per fold!
        fold_model = clone(ml_model)
        # Training

        # Convert the fingerprint and the label to a list
        train_x = df.iloc[train_index].fp.tolist()
        train_y = df.iloc[train_index].active.tolist()

        # Fit the model
        fold_model.fit(train_x, train_y)

        # Testing

        # Convert the fingerprint and the label to a list
        test_x = df.iloc[test_index].fp.tolist()
        test_y = df.iloc[test_index].active.tolist()

        # Performance for each fold
        accuracy, sens, spec, auc = model_performance(fold_model, test_x, test_y, verbose)

        # Save results
        acc_per_fold.append(accuracy)
        sens_per_fold.append(sens)
        spec_per_fold.append(spec)
        auc_per_fold.append(auc)

    # Print statistics of results
    print(
        f"Mean accuracy: {np.mean(acc_per_fold):.2f} \t"
        f"and std : {np.std(acc_per_fold):.2f} \n"
        f"Mean sensitivity: {np.mean(sens_per_fold):.2f} \t"
        f"and std : {np.std(sens_per_fold):.2f} \n"
        f"Mean specificity: {np.mean(spec_per_fold):.2f} \t"
        f"and std : {np.std(spec_per_fold):.2f} \n"
        f"Mean AUC: {np.mean(auc_per_fold):.2f} \t"
        f"and std : {np.std(auc_per_fold):.2f} \n"
        f"Time taken : {time.time() - t0:.2f}s\n"
    )

    return acc_per_fold, sens_per_fold, spec_per_fold, auc_per_fold<jupyter_output><empty_output><jupyter_text>**Cross-validation**We now apply cross-validation and show the statistics for all three ML models. In real world conditions, cross-validation usually applies 5 or more folds, but for the sake of performance we will reduce it to 3. You can change the value of `N_FOLDS` in this cell below.<jupyter_code>N_FOLDS = 3<jupyter_output><empty_output><jupyter_text>_Note_: Next cell takes long to execute<jupyter_code>for model in models:
    print("\n======= ")
    print(f"{model['label']}")
    crossvalidation(model["model"], compound_df, n_folds=N_FOLDS)<jupyter_output>======= 
Model_RF
Mean accuracy: 0.82 	and std : 0.00 
Mean sensitivity: 0.89 	and std : 0.01 
Mean specificity: 0.70 	and std : 0.01 
Mean AUC: 0.88 	and std : 0.01 
Time taken : 0.70s


======= 
Model_SVM
Mean accuracy: 0.82 	and std : 0.00 
Mean sensitivity: 0.90 	and std : 0.00 
Mean specificity: 0.69 	and std : 0.02 
Mean AUC: 0.89 	and std : 0.01 
Time taken : 5.16s


======= 
Model_ANN
Mean accuracy: 0.79 	and std : 0.01 
Mean sensitivity: 0.88 	and std : 0.00 
Mean specificity: 0.65 	and std : 0.02 
Mean AUC: 0.85 	and std : 0.01 
Time taken : 2.15s<jupyter_text>We look at the cross-validation performance for molecules encoded using Morgan fingerprint and not MACCS keys.<jupyter_code># Reset data frame
compound_df = chembl_df.copy()
# Use Morgan fingerprint with radius 3
compound_df["fp"] = compound_df["smiles"].apply(smiles_to_fp, args=("morgan3",))
compound_df.head(3)
# NBVAL_CHECK_OUTPUT<jupyter_output><empty_output><jupyter_text>_Note_: Next cell takes long to execute<jupyter_code>for model in models:
    if model["label"] == "Model_SVM":
        # SVM is super slow with long fingerprints
        # and will have a performance similar to RF
        # We can skip it in this test, but if you want
        # to run it, feel free to replace `continue` with `pass`
        continue
    print("\n=======")
    print(model["label"])
    crossvalidation(model["model"], compound_df, n_folds=N_FOLDS)<jupyter_output>=======
Model_RF
Mean accuracy: 0.83 	and std : 0.00 
Mean sensitivity: 0.91 	and std : 0.01 
Mean specificity: 0.71 	and std : 0.02 
Mean AUC: 0.90 	and std : 0.01 
Time taken : 2.29s


=======
Model_ANN
Mean accuracy: 0.81 	and std : 0.01 
Mean sensitivity: 0.87 	and std : 0.01 
Mean specificity: 0.70 	and std : 0.01 
Mean AUC: 0.87 	and std : 0.01 
Time taken : 11.42s<jupyter_text>Regresion modelsSo far, we have been running classification models (is the compound going to be active or not, indicated by the active column we added in the dataframe). The next step is to see if we can also train a regression model, which can be particularly interesting if we want to design our own ligands (e.g. do we predict the compound to actually bind better). We will remove the censored data (the ones with no pChEMBL value), and train another model. Make a new function that works with regression:<jupyter_code>def crossvalidation_reg(ml_model, df, n_folds=5, verbose=False):
    """
    Machine learning model training and validation in a cross-validation loop.

    Parameters
    ----------
    ml_model: sklearn model object
        The machine learning model to train.
    df: pd.DataFrame
        Data set with SMILES and their associated activity labels.
    n_folds: int, optional
        Number of folds for cross-validation.
    verbose: bool, optional
        Performance measures are printed.

    Returns
    -------
    None

    """
    t0 = time.time()
    # Shuffle the indices for the k-fold cross-validation
    kf = KFold(n_splits=n_folds, shuffle=True, random_state=SEED)

    # Results for each of the cross-validation folds
    MAE_per_fold = []
    RMSE_per_fold = []

    # Loop over the folds
    for train_index, test_index in kf.split(df):
        # clone model -- we want a fresh copy per fold!
        fold_model = clone(ml_model)
        # Training

        # Convert the fingerprint and the label to a list
        train_x = df.iloc[train_index].fp.tolist()
        train_y = df.iloc[train_index].pKi.tolist()

        # Fit the model
        fold_model.fit(train_x, train_y)

        # Testing

        # Convert the fingerprint and the label to a list
        test_x = df.iloc[test_index].fp.tolist()
        test_y = df.iloc[test_index].pKi.tolist()
        
        test_results = fold_model.predict(test_x)
        # Prediction probability on test set
        from sklearn import metrics

        MAE_per_fold.append(metrics.mean_absolute_error(test_y, test_results))
        #print('Mean Squared Error (MSE):', metrics.mean_squared_error(test_y, test_results))
        RMSE_per_fold.append(np.sqrt(metrics.mean_squared_error(test_y, test_results)))
        #mape = np.mean(np.abs((gt - pred) / np.abs(gt)))
        #print('Mean Absolute Percentage Error (MAPE):', round(mape * 100, 2))
        #print('Accuracy:', round(100*(1 - mape), 2))
    return(MAE_per_fold,RMSE_per_fold,fold_model)<jupyter_output><empty_output><jupyter_text>Next, create a new dataset (without censored data) and run.<jupyter_code>fingerprint_to_model = compound_df.fp.tolist()
label_to_model = compound_df.pKi.tolist()

# Split data randomly in train and test set
# note that we use test/train_x for the respective fingerprint splits
# and test/train_y for the respective label splits
(
    static_train_x,
    static_test_x,
    static_train_y,
    static_test_y,
) = train_test_split(fingerprint_to_model, label_to_model, test_size=0.2, random_state=12)
splits = [static_train_x, static_test_x, static_train_y, static_test_y]
# NBVAL_CHECK_OUTPUT
print("Training data size:", len(static_train_x))
print("Test data size:", len(static_test_x))
# Kick out NaN in the pChEMBL values
compound_df_noNaN = compound_df.dropna()

# Append RF Regressor model
#models.append({"label": "Model_RF_reg", "model": RandomForestRegressor})
#Train model with RandomForestRegressor
regressor = RandomForestRegressor()
MAE, RMSE,trained_model = crossvalidation_reg(regressor , compound_df_noNaN, n_folds=N_FOLDS)

print(
f"Mean Absolute Error (MAE): {np.mean(MAE):.2f} \t"
f"and std : {np.std(MAE):.2f} \n"
f"Root Mean Square Error (RMSE): {np.mean(RMSE):.2f} \t"
f"and std : {np.std(RMSE):.2f} \n"
)<jupyter_output>Mean Absolute Error (MAE): 0.51 	and std : 0.02 
Root Mean Square Error (RMSE): 0.69 	and std : 0.03
{'fraction_non_alphanumeric': 0.09321436957032167, 'fraction_numerical': 0.018247073424345086, 'mean_word_length': 4.128849131257526, 'pattern_counts': {'":': 14, '<': 84, '<?xml version=': 0, '>': 86, 'https://': 21, 'lorem ipsum': 0, 'www.': 2, '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': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '423063', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 9677, 'n_tokens_neox': 8800, 'n_words': 3083, 'path': 'teachopencadd/talktorials/T007_compound_activity_machine_learning/talktorial.ipynb', 'repo_name': 'jesperswillem/teachopencadd_a2a', 'size': '29894'}