LinearRegressionMultipleModels4CDK2_Ki.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "# **Multiple Linear Regression Models for the CDK2 K<sub>i</sub> Dataset**\n",
        "\n",
        "This code builds multiple regression models using Ordinary Least Squares ([Géron, 2023](https://www.isbns.net/isbn/9781098125974/); [Nield, 2022](https://www.isbns.net/isbn/9781098102937/)) for data from the CDK2 K<sub>i</sub> dataset ([de Azevedo et al., 2024](https://doi.org/10.1002/jcc.27449)). It employs the [Scikit-Learn](https://scikit-learn.org/stable/) library ([Pedregosa et al., 2011](https://doi.org/10.48550/arXiv.1201.0490)) to generate the regression model. It evaluates the predictive performance of the regression model using metrics recommended by [Walsh et al., 2021](https://doi.org/10.1038/s41592-021-01205-4) and de [Azevedo et al., 2024](https://doi.org/10.1002/jcc.27449).\n",
        "\n",
        "<br></br>\n",
        "<img src=\"https://drive.usercontent.google.com/download?id=1URfpJAcHFpwsqCIQkbCvpmArX-t34Sua&export=view&authuser=0\" width=320 alt=\"STS-95\">\n",
        "<br></br>\n",
        "Protein-ligand structure used to train a model to predict binding affinity.\n",
        "<br></br>\n",
        "**References**\n",
        "\n",
        "de Azevedo WF Jr, Quiroga R, Villarreal MA, da Silveira NJF,\n",
        "Bitencourt-Ferreira G, da Silva AD, Veit-Acosta M, Oliveira PR, Tutone M, Biziukova N, Poroikov V, Tarasova O, Baud S. SAnDReS 2.0: Development of machine-learning models to explore the scoring function space. J Comput Chem. 2024; 45(27): 2333-2346.\n",
        "PMID: 38900052. [DOI: 10.1002/jcc.27449](https://doi.org/10.1002/jcc.27449) [PubMed](https://pubmed.ncbi.nlm.nih.gov/38900052/)\n",
        "<br></br>\n",
        "Géron, A. Hands-on Machine Learning with Scikit-Learn, Keras, and TensorFlow 3e: Concepts, Tools, and Techniques to Build Intelligent Systems, 3rd ed.; O’Reilly Media: Sebastopol, CA, 2023.\n",
        "[ISBN: 978-1-098-12597-4](https://www.isbns.net/isbn/9781098125974/)\n",
        "<br></br>\n",
        "Nield, T. Essential Math for Data Science: Take Control of Your Data with Fundamental Linear Algebra, Probability, and Statistics; O’Reilly Media: Sebastopol, CA, 2022. [ISBN: 978-1-098-10293-7](https://www.isbns.net/isbn/9781098102937/)\n",
        "<br></br>\n",
        "Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, Blondel M, Prettenhofer P, Weiss R, Dubourg V, Verplas J, Passos A, Cournapeau D, Brucher M, Perrot M, Duchesnay E. Scikitlearn: Machine Learning in Python. J Mach Learn Res., 2011; 12:2825–2830. [DOI: 10.48550/arXiv.1201.0490](https://doi.org/10.48550/arXiv.1201.0490)\n",
        "<br></br>\n",
        "Walsh I, Fishman D, Garcia-Gasulla D, Titma T, Pollastri G; ELIXIR Machine Learning Focus Group; Harrow J, Psomopoulos FE, Tosatto SCE. DOME: recommendations for supervised machine learning validation in biology. Nat Methods. 2021; 18(10): 1122-1127. [DOI: 10.1038/s41592-021-01205-4](https://doi.org/10.1038/s41592-021-01205-4) [PubMed](https://pubmed.ncbi.nlm.nih.gov/34316068/)\n",
        "<br></br>\n",
        "\n",
        "It follows the code.\n"
      ],
      "metadata": {
        "id": "LbZUSKh33KMF"
      },
      "id": "LbZUSKh33KMF"
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "83f24dc7-816d-4be2-8b39-4d0da983b3b4",
      "metadata": {
        "id": "83f24dc7-816d-4be2-8b39-4d0da983b3b4",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "070e2576-314d-4116-a2d7-bf1f5fa34205"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "\n",
            "Downloading and preprocessing file: cdk2_ki_training_set.csv...done!\n",
            "\n",
            "Downloading and preprocessing file: cdk2_ki_test_set.csv...done!\n",
            "\n",
            "Generating regression models...done!\n",
            "\n",
            "Number of generated regression models:  924\n",
            "\n",
            "Metrics saved to file: regression_metrics.csv\n"
          ]
        }
      ],
      "source": [
        "#!/usr/bin/env python3\n",
        "#\n",
        "################################################################################\n",
        "# Dr. Walter F. de Azevedo, Jr.                                                #\n",
        "# [Scopus](https://www.scopus.com/authid/detail.uri?authorId=7006435557)       #\n",
        "# [GitHub](https://github.com/azevedolab)                                      #\n",
        "# July 20, 2024                                                                #\n",
        "################################################################################\n",
        "#\n",
        "################################################################################\n",
        "# Import section                                                               #\n",
        "################################################################################\n",
        "import pandas as pd\n",
        "from itertools import combinations\n",
        "from matplotlib import pyplot as plt\n",
        "from sklearn import preprocessing\n",
        "from sklearn.linear_model import LinearRegression\n",
        "\n",
        "################################################################################\n",
        "# Define cross_validation() function                                           #\n",
        "################################################################################\n",
        "# Function to set up k-fold class. Kfold class to build a n fold\n",
        "# cross-validation loop and test the generalization ability of regression. With\n",
        "# cross-validation, we generally obtain a more conservative estimate(that is,\n",
        "# the error is larger). The cross-validation estimate is a better estimate of\n",
        "# how well we could generalize to predict on unseen data.\n",
        "#\n",
        "################################################################################\n",
        "# Reference                                                                    #\n",
        "################################################################################\n",
        "# Coelho LP, Richert W. (2015) Building Machine Learning Systems with\n",
        "# Python. 2nd ed. Packt Publishing Ltd. Birmingham UK. 301 pp. See page 162\n",
        "# (Cross-validation for regression)\n",
        "def cross_validation(model,X,y,n_splits,random_state,verbose):\n",
        "\n",
        "    # Import section\n",
        "    from sklearn.model_selection import KFold, cross_val_score\n",
        "    from warnings import simplefilter\n",
        "    import numpy as np\n",
        "\n",
        "    # Set up k-fold class\n",
        "    kf = KFold(n_splits=n_splits,shuffle=True, random_state=random_state)\n",
        "\n",
        "    # Ignore all future warnings\n",
        "    simplefilter(action='ignore', category=DeprecationWarning)\n",
        "\n",
        "    # Looping through kf.split()\n",
        "    for train,test in kf.split(X):\n",
        "\n",
        "        # Generate regression model\n",
        "        model.fit(X[train],y[train])\n",
        "\n",
        "    # Show Walsh metrics if requestes\n",
        "    if verbose:\n",
        "        # Show average coefficient of determination using n-fold crossvalidation\n",
        "        scores = cross_val_score(model,X,y,cv=kf)\n",
        "        msg_out = \"Average coefficient of determination using n-fold \"\n",
        "        msg_out += \"cross-validation\"\n",
        "        print(\"\\n\"+msg_out+\":\",np.mean(scores))\n",
        "\n",
        "    # Return model\n",
        "    return model\n",
        "\n",
        "################################################################################\n",
        "# Define Walsh_and_de_Azevedo_regression_metrics() function                    #\n",
        "################################################################################\n",
        "def Walsh_and_de_Azevedo_regression_metrics(y,y_pred):\n",
        "    \"\"\"Function to calculate regression metrics as recommended by\n",
        "    [Walsh et 2021](https://pubmed.ncbi.nlm.nih.gov/34316068/)\n",
        "    and\n",
        "    [de Azevedo Jr et al., 2024](https://pubmed.ncbi.nlm.nih.gov/38900052/))\n",
        "\n",
        "    Parameters:\n",
        "    y: array with experimental data (dimensions: (n,) )\n",
        "    y_pred: array with predicted values (dimensions: (n,) )\n",
        "\n",
        "    where n is the number of instances\n",
        "    \"\"\"\n",
        "\n",
        "    # Import section for this function\n",
        "    from sklearn.metrics import median_absolute_error\n",
        "    from scipy import stats\n",
        "    import numpy as np\n",
        "\n",
        "    # Determine metrics\n",
        "    n = len(y)\n",
        "    y_bar = y.mean()\n",
        "    rss = np.sum( np.square(y_pred - y) )\n",
        "    rmse =  np.sqrt( (np.sum( (y - y_pred)**2 ))/n )\n",
        "    ss_res = np.sum( (y - y_pred)**2 )\n",
        "    ss_tot = np.sum( (y - y_bar)**2 )\n",
        "    R2 = 1 - ss_res/ss_tot\n",
        "    mae = median_absolute_error(y,y_pred)\n",
        "    r_pearson,p_pearson = stats.pearsonr(y,y_pred)\n",
        "    rho,p_spearman = stats.spearmanr(y,y_pred)\n",
        "    core_dome = rmse**2 + mae**2 + (R2 -1)**2\n",
        "    dome =    np.sqrt( core_dome )\n",
        "    edomer2 = np.sqrt( core_dome + (r_pearson**2 - 1)**2 )\n",
        "    edomerho =np.sqrt( core_dome + (rho**2 - 1)**2 )\n",
        "    edome =   np.sqrt( core_dome + (r_pearson**2 - 1)**2 + (rho**2 - 1)**2 )\n",
        "\n",
        "    # Return metrics\n",
        "    return r_pearson,p_pearson,rho,p_spearman,rss,rmse,mae,R2,dome,edomer2,\\\n",
        "    edomerho,edome\n",
        "\n",
        "################################################################################\n",
        "# Training set                                                                 #\n",
        "################################################################################\n",
        "file_in = \"cdk2_ki_training_set.csv\"             # Update here for a new dataset\n",
        "print(\"\\nDownloading and preprocessing file: \"+file_in,end=\"...\")\n",
        "url = \"https://bit.ly/CDK2_Ki_Training_Set\"      # Update here for a new dataset\n",
        "feature_list =[\"Gauss 2\",\"C\",\"Q\",\"Torsions\",\"O\"] # Update here for a new dataset\n",
        "training_set_data = pd.read_csv(url,delimiter=\",\")\n",
        "# X_train = training_set_data[feature_list]\n",
        "y_train = training_set_data.iloc[:,8]            # Update here for a new dataset\n",
        "# scaler_train = preprocessing.StandardScaler().fit(X_train)\n",
        "# X_train = scaler_train.transform(X_train)\n",
        "print(\"done!\")\n",
        "\n",
        "################################################################################\n",
        "# Test set                                                                     #\n",
        "################################################################################\n",
        "file_in = \"cdk2_ki_test_set.csv\"                 # Update here for a new dataset\n",
        "print(\"\\nDownloading and preprocessing file: \"+file_in,end=\"...\")\n",
        "url = \"https://bit.ly/CDK2_Ki_Test_Set\"          # Update here for a new dataset\n",
        "test_set_data = pd.read_csv(url,delimiter=\",\")\n",
        "# X_test = test_set_data[feature_list]             # Update here for a new dataset\n",
        "y_test = test_set_data.iloc[:,8]                 # Update here for a new dataset\n",
        "# scaler_test = preprocessing.StandardScaler().fit(X_test)\n",
        "# X_test = scaler_test.transform(X_test)\n",
        "print(\"done!\")\n",
        "\n",
        "################################################################################\n",
        "# Define additional inputs                                                     #\n",
        "################################################################################\n",
        "plt_scatter_X_label_str,plt_scatter_y_label_str = \"pK$_i$ (Experimental)\",\\\n",
        "                        \"pK$_i$ (Predicted)\"            # Define labels for axes\n",
        "plt_dpi = 1500                                          # Plot dpi\n",
        "\n",
        "################################################################################\n",
        "# Build regression models with Ordinary Least Squares (LinearRegression)       #\n",
        "################################################################################\n",
        "# Set up initial values\n",
        "n_models = 0\n",
        "rmse_list = []\n",
        "l = 6\n",
        "y_list = [y_train,y_test]\n",
        "data_list = [\"Training Set\",\"Test Set\"]\n",
        "lst_in = [\"Gauss 1\",\"Gauss 2\",\"Repulsion\",\"Hydrophobic\",\"Hydrogen\",\"Torsional\",\n",
        "\"Torsions\",\"Q\",\"Average Q\",\"C\",\"N\",\"O\"]\n",
        "data_out = \"Features,Data,n,r,p-value(r),r2,rho,p-value(rho),RSS,RMSE,MAE,R2,\"\n",
        "data_out += \"DOME,EDOMEr2,EDOMErho,EDOME\\n\"\n",
        "\n",
        "# Get all combinations of elements in lst_in and length l\n",
        "comb_lst_in = combinations(lst_in, l)\n",
        "\n",
        "# Show message\n",
        "msg_out = \"\\nGenerating regression models\"\n",
        "print(msg_out,end=\"...\")\n",
        "\n",
        "# Looping through comb_lst_in\n",
        "for feature_list in list(comb_lst_in):\n",
        "\n",
        "  # Make it a list\n",
        "  feature_list = list(feature_list)\n",
        "\n",
        "  try:\n",
        "    # For training set\n",
        "    X_train = training_set_data[feature_list]\n",
        "    scaler_train = preprocessing.StandardScaler().fit(X_train)\n",
        "    X_train = scaler_train.transform(X_train)\n",
        "\n",
        "    # For test set\n",
        "    X_test = test_set_data[feature_list]\n",
        "    scaler_test = preprocessing.StandardScaler().fit(X_test)\n",
        "    X_test = scaler_test.transform(X_test)\n",
        "    n_models += 1\n",
        "\n",
        "    # Instantiate an object of LinearRegression() class\n",
        "    lin_reg = LinearRegression().fit(X_train, y_train)\n",
        "\n",
        "    # Call cross_validation() function\n",
        "    lin_reg = cross_validation(lin_reg,X_train,y_train,5,46,False)\n",
        "    y_train_pred_reg = lin_reg.predict(X_train)\n",
        "    y_test_pred_reg = lin_reg.predict(X_test)\n",
        "    y_pred_list = [y_train_pred_reg,y_test_pred_reg]\n",
        "\n",
        "    # Looping through data\n",
        "    for count,y in enumerate(y_list):\n",
        "      # Call Walsh_and_de_Azevedo_regression_metrics() function\n",
        "      r_pearson,p_pearson,rho,p_spearman,rss,rmse,mae,R2,dome,edomer2,\\\n",
        "      edomerho,edome = \\\n",
        "      Walsh_and_de_Azevedo_regression_metrics(y_list[count],\n",
        "                                                            y_pred_list[count])\n",
        "\n",
        "      # Set up output line\n",
        "      line_out = str(feature_list).replace(\"'\",\"\").replace(\"[\",\"\").\\\n",
        "                                                replace(\"]\",\"\").replace(\",\",\" \")\n",
        "      line_out += \",\"+data_list[count]+\",\"+str(len(y))+\",\"\n",
        "      line_out += \"{:.4f}\".format(r_pearson)+\",\"+\"{:.4e}\".format(p_pearson)+\",\"\n",
        "      line_out += \"{:.4f}\".format(r_pearson**2)+\",\"+\"{:.4f}\".format(rho)+\",\"\n",
        "      line_out += \"{:.4e}\".format(p_spearman)+\",\"+\"{:.4e}\".format(rss)+\",\"\n",
        "      line_out += \"{:.4f}\".format(rmse)+\",\"+\"{:.4f}\".format(mae)+\",\"\n",
        "      line_out += \"{:.4f}\".format(R2)+\",\"+\"{:.4f}\".format(dome)+\",\"\n",
        "      line_out += \"{:.4f}\".format(edomer2)+\",\"+\"{:.4f}\".format(edomerho)+\",\"\n",
        "      line_out += \"{:.4f}\".format(edome)\n",
        "\n",
        "      # Update data_out\n",
        "      data_out += line_out+\"\\n\"\n",
        "\n",
        "  except:\n",
        "    print(\"\\nError I didn't generate a regression model for \",feature_list)\n",
        "\n",
        "print(\"done!\")\n",
        "\n",
        "# Open a new file and write data_out\n",
        "fo_out = open(\"regression_metrics.csv\",\"w\")\n",
        "fo_out.write(data_out)\n",
        "fo_out.close()\n",
        "print(\"\\nNumber of generated regression models: \",n_models)\n",
        "print(\"\\nMetrics saved to file: regression_metrics.csv\")\n",
        "################################################################################"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.12.4"
    },
    "colab": {
      "provenance": []
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}