CNN_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### This notebook performs regression using a Convolutional Neural Net on molecular descriptors.\n",
    "#### There are 1119 molecular descriptors for each molecule. A dummy zero descriptor is added to make the total number of descriptors 1120.\n",
    "#### Then the descriptors are reshaped into 35 by 32 (35*32)  grayscale images.\n",
    "#### Using Leave one out cross-validation data generated by create_dataset.ipynb, models are trained on different folds of training data to regress on pIC50 values.\n",
    "#### The trained models are saved in models_info folder.\n",
    "#### Various parameters like R2, tropsha's criteria etc are also evaluated for each model and saved in models_info/logs/logs.txt\n",
    "#### Please note that the above mentioned folders are empty. They will be populated once this notebook is run.\n",
    "#### Best models based on R2 test are saved in best_models folder. (.json files are model architecture and .h5 files re model weights) Performance parameters for best chosen models are available in best_models/performance_parameters_of_best_models/performance_parameters.xlsx. \n",
    "#### They were evaluated through XternalValidationPlus (https://sites.google.com/site/dtclabxvplus/) by passing the test set predictions of each model. These predictions are available in models_info/data_to_evaluate_xternal_validation_parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import sys\n",
    "from keras.models import Sequential,model_from_json\n",
    "from keras.layers import Dense, Conv2D, Flatten, MaxPooling2D,UpSampling2D, Dropout\n",
    "from keras.callbacks import EarlyStopping,ModelCheckpoint\n",
    "pd.set_option('display.max_rows', None)\n",
    "import matplotlib.pyplot as plt\n",
    "import tensorflow as tf\n",
    "from sklearn.metrics import r2_score\n",
    "from sklearn.metrics import mean_absolute_error\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.linear_model import LinearRegression\n",
    "import math\n",
    "import pickle\n",
    "import logging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Evaluate tropsha's criteria\n",
    "def tropsha(predicted_holdout_values,holdout_labels):\n",
    "    \n",
    "    \n",
    "    # Find attributes of test set\n",
    "    reg = LinearRegression().fit(np.reshape(predicted_holdout_values,(-1,1)),holdout_labels)\n",
    "    \n",
    "\n",
    "    r2 = reg.score(np.reshape(predicted_holdout_values,(-1,1)),holdout_labels)\n",
    "    \n",
    "    reg = LinearRegression(fit_intercept = False).fit(np.reshape(predicted_holdout_values,(-1,1)),holdout_labels)\n",
    "\n",
    "    r02 = reg.score(np.reshape(predicted_holdout_values,(-1,1)),holdout_labels)\n",
    "    \n",
    "    k = reg.coef_\n",
    "\n",
    "    r0_dash2 = reg.score(np.reshape(predicted_holdout_values,(-1,1)),np.reshape(np.array(holdout_labels),(-1,1)))\n",
    "\n",
    "    k_dash = reg.coef_\n",
    "    \n",
    "    print(\"R2: (>0.6): \"+str(r2))\n",
    "    print(\"R02: \"+str(r02))\n",
    "    print(\"R0_dash2: \"+str(r0_dash2))\n",
    "    print(\"R02-R0_dash2: (<0.3): \"+str(np.abs(r02 - r0_dash2)))\n",
    "    print(\"k: (0.85,1.15): \"+str(k))\n",
    "    print(\"(R2-R02)/R2: (<0.1)\"+str((r2-r02)/r2))\n",
    "    print(\"k_dash: (0.85,1.15): \"+str(k_dash))\n",
    "    print(\"(R2-R0_dash2)/R2 (<0.1): \"+str((r2-r0_dash2)/r2))\n",
    "    if (r2 > 0.6 and np.abs(r02 - r0_dash2)<0.3 and k>=0.85 and k<=1.15 and (r2-r02)/r2<0.1 and k_dash>=0.85 and k_dash<=1.15\n",
    "        and (r2-r0_dash2)/r2<0.1):\n",
    "        return True\n",
    "    else:\n",
    "        return False\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read test set\n",
    "csv_df_test = pd.read_csv('data/test_compounds.csv')\n",
    "csv_df_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# output folder for saving models and weights\n",
    "out_folder_name = \"models_info/\"\n",
    "\n",
    "# Log location for logging model performance\n",
    "logging.basicConfig(filename=out_folder_name+\"logs/log.txt\", level=logging.INFO)\n",
    "\n",
    "# No of models due to data folds generated by Leave-one-out cross validation. In this notebook there are 69 datapoints \n",
    "# in the training set. Hence n_models = 69\n",
    "n_models = 94-len(csv_df_test)\n",
    "# Image dimensions for molecular descriptors reshaping\n",
    "img_height = 35\n",
    "img_width = 32\n",
    "\n",
    "# Load scaler used to reshape training data. This will be used to inverse transform data to get actual pIC50 values. \n",
    "# Note that the data used to train CNN is already min-max scaled between 0 and 1.\n",
    "scaler = pickle.load(open('scaler_data/scaler.dat','rb'))\n",
    "test_data = csv_df_test.loc[:,'nAcid':'values']\n",
    "\n",
    "# min max inverse transform for test data\n",
    "test_data_inverse_transform_actual = scaler.inverse_transform(test_data)\n",
    "test_values_inverse_transform_actual = test_data_inverse_transform_actual[:,len(test_data_inverse_transform_actual[0])-1]\n",
    "test_data = test_data.loc[:,'nAcid':'Zagreb']\n",
    "test_labels = np.reshape(np.array(csv_df_test['values']),-1)\n",
    "num_of_test_examples = len(test_labels)\n",
    "\n",
    "# Add a dummy 0 to increase descriptor columns to 1120 so that descriptors can be reshaped into 35 by 32 grayscale images.\n",
    "test_data['zero'] = 0\n",
    "test_imgs = np.reshape(np.array(test_data),(num_of_test_examples,img_height,img_width,1))\n",
    "\n",
    "r2_training_list = []\n",
    "r2_validation_list = []\n",
    "dict_of_metrics = {}\n",
    "\n",
    "# looping and creating models\n",
    "for i in range(0,n_models):\n",
    "    # Read training data\n",
    "    csv_df_training = pd.read_csv('data/training_set_'+str(i)+'.csv')\n",
    "    # Read leave one out data\n",
    "    csv_df_loo = pd.read_csv('data/loo_set__'+str(i)+'.csv')\n",
    "\n",
    "    training_data = csv_df_training.loc[:,'nAcid':'values']\n",
    "    loo_data = csv_df_loo.loc[:,'nAcid':'values']\n",
    "    num_of_training_examples = len(training_data)\n",
    "    num_of_loo_examples = len(loo_data)\n",
    "    \n",
    "    # Minmax scaler inverse transform for  train/loo set\n",
    "    training_data_inverse_transform_actual = scaler.inverse_transform(training_data)\n",
    "    loo_data_inverse_transform_actual = scaler.inverse_transform(loo_data)\n",
    "    \n",
    "    training_values_inverse_transform_actual = training_data_inverse_transform_actual[:,len(training_data_inverse_transform_actual[0])-1]\n",
    "    loo_values_inverse_transform_actual = loo_data_inverse_transform_actual[:,len(loo_data_inverse_transform_actual[0])-1]\n",
    "    ######\n",
    "    training_labels = np.reshape(np.array(training_data['values']),-1)\n",
    "    loo_labels = loo_data['values']\n",
    "    \n",
    "    training_data = training_data.loc[:,'nAcid':'Zagreb']\n",
    "    loo_data = loo_data.loc[:,'nAcid':'Zagreb']\n",
    "    \n",
    "    # Add a dummy 0 to increase descriptor columns to 1120 so that descriptors can be reshaped into 35 by 32 grayscale images.\n",
    "    training_data['zero'] = 0\n",
    "    loo_data['zero'] = 0\n",
    "    \n",
    "    training_imgs = np.reshape(np.array(training_data),(num_of_training_examples,img_height,img_width,1))\n",
    "    loo_imgs = np.reshape(np.array(loo_data),(num_of_loo_examples,img_height,img_width,1))\n",
    "    \n",
    "    # CNN\n",
    "    #create model\n",
    "    model = Sequential()\n",
    "    #add model layers\n",
    "    model.add(Conv2D(4, kernel_size=5, activation='relu', input_shape=(img_height,img_width,1),padding='valid'))\n",
    "    model.add(Dropout(0.02))\n",
    "    model.add(Conv2D(6, kernel_size=5, activation='relu',padding='valid'))\n",
    "    model.add(Dropout(0.04))\n",
    "    model.add(Conv2D(8, kernel_size=5, activation='relu',padding='valid'))\n",
    "    model.add(Dropout(0.06))\n",
    "    model.add(Flatten())\n",
    "    model.add(Dense(128, activation='relu'))\n",
    "    model.add(Dropout(0.2))\n",
    "    model.add(Dense(128, activation='relu'))\n",
    "    model.add(Dropout(0.2))\n",
    "    model.add(Dense(128, activation='relu'))\n",
    "    model.add(Dropout(0.2))\n",
    "    model.add(Dense(128, activation='relu'))\n",
    "    model.add(Dropout(0.1))\n",
    "    model.add(Dense(1, activation='relu'))\n",
    "    model.compile(optimizer='adam', loss='mse')\n",
    "    \n",
    "    # Check point for saving best model\n",
    "    mc = ModelCheckpoint(out_folder_name+\"/models/model_\"+str(i)+\".h5\", monitor='val_loss', mode='min', verbose=0, save_best_only=True,\n",
    "                         save_weights_only=True)\n",
    "      \n",
    "    \n",
    "    print(\"Model-\"+str(i))\n",
    "    logging.info(\"Model-\"+str(i))\n",
    "    \n",
    "    # train the model\n",
    "    history = model.fit(training_imgs, training_labels,\n",
    "                        validation_data=(test_imgs,test_labels), \n",
    "                        epochs=2000,batch_size=25,verbose = 0,callbacks=[mc])\n",
    "    training_loss_list = history.history['loss']\n",
    "    validation_loss_list = history.history['val_loss']\n",
    "    \n",
    "     # serialize model to JSON\n",
    "    model_json = model.to_json()\n",
    "    with open(out_folder_name+\"/models/model_\"+str(i)+\".json\", \"w\") as json_file:\n",
    "        json_file.write(model_json)\n",
    "    \n",
    "    # Load  model\n",
    "    json_file = open(out_folder_name+'/models/model_'+str(i)+'.json')\n",
    "    loaded_model_json = json_file.read()\n",
    "    json_file.close()\n",
    "    loaded_model = model_from_json(loaded_model_json)\n",
    "    loaded_model.load_weights(out_folder_name+'/models/model'+str(i)+'.h5')\n",
    "    \n",
    "    # Predict training set\n",
    "    predicted_train_values = loaded_model.predict(training_imgs)\n",
    "    # min max scaler inverse transform for predicted train values\n",
    "    training_data_predicted = training_data.loc[:,'nAcid':'Zagreb']\n",
    "    training_data_predicted['values'] = predicted_train_values\n",
    "    training_data_inverse_transform_predicted = scaler.inverse_transform(training_data_predicted)\n",
    "    training_value_inverse_transform_predicted = training_data_inverse_transform_predicted[:,len(training_data_inverse_transform_predicted[0])-1]\n",
    "    \n",
    "    # Predict loo set\n",
    "    predicted_loo_values = loaded_model.predict(loo_imgs)\n",
    "    # min max scaler inverse transform for predicted loo values\n",
    "    loo_data_predicted = loo_data.loc[:,'nAcid':'Zagreb']\n",
    "    loo_data_predicted['values'] = predicted_loo_values\n",
    "    loo_data_inverse_transform_predicted = scaler.inverse_transform(loo_data_predicted)\n",
    "    loo_value_inverse_transform_predicted = loo_data_inverse_transform_predicted[:,len(loo_data_inverse_transform_predicted[0])-1]\n",
    "    \n",
    "    # Predict test set\n",
    "    predicted_test_values = loaded_model.predict(test_imgs)\n",
    "    test_data_predicted = test_data.loc[:,'nAcid':'Zagreb']\n",
    "    test_data_predicted['values'] = predicted_test_values\n",
    "    # min max scaler inverse transform for predicted test values\n",
    "    test_data_inverse_transform_predicted = scaler.inverse_transform(test_data_predicted)\n",
    "    test_value_inverse_transform_predicted = test_data_inverse_transform_predicted[:,len(test_data_inverse_transform_predicted[0])-1]\n",
    "    \n",
    "    predicted_test_values = np.reshape(predicted_test_values,-1)\n",
    "    \n",
    "    # Calculate model r2\n",
    "    r2_training = r2_score(training_values_inverse_transform_actual,np.reshape(training_value_inverse_transform_predicted,-1))\n",
    "    r2_test = r2_score(test_values_inverse_transform_actual,test_value_inverse_transform_predicted)\n",
    "    \n",
    "    # Calculate LOO-Q2\n",
    "    avg_predicted_train_values = np.mean(training_value_inverse_transform_predicted)\n",
    "    avg_train_values = np.mean(np.array(training_values_inverse_transform_actual))\n",
    "    TSS = np.mean(np.square(np.reshape(training_values_inverse_transform_actual,-1)-avg_train_values))\n",
    "    PRESS = np.square(np.reshape(loo_value_inverse_transform_predicted,-1)-np.reshape(loo_values_inverse_transform_actual,-1))\n",
    "    LOO_Q2 = 1-(PRESS/TSS)\n",
    "    \n",
    "    # Calculate Q2-ext-F1\n",
    "    q2_ext_numerator = np.sum(np.square(test_value_inverse_transform_predicted - np.reshape(np.array(test_values_inverse_transform_actual),-1)))\n",
    "    q2_ext_denominator = np.sum(np.square(np.reshape(np.array(test_values_inverse_transform_actual),-1) - avg_train_values))\n",
    "    q2_ext = 1-(q2_ext_numerator/q2_ext_denominator)\n",
    "     \n",
    "    mae_test = mean_absolute_error(test_values_inverse_transform_actual, test_value_inverse_transform_predicted)\n",
    "    se_training = math.sqrt(mean_squared_error(training_values_inverse_transform_actual,training_value_inverse_transform_predicted))\n",
    "    mae_training = mean_absolute_error(training_values_inverse_transform_actual, training_value_inverse_transform_predicted)\n",
    "    max_training_response_value = np.max(training_values_inverse_transform_actual)\n",
    "    min_training_response_value = np.min(training_values_inverse_transform_actual)\n",
    "    \n",
    "    print(\"Training response avg: \"+str(avg_train_values))\n",
    "    logging.info(\"Training response avg: \"+str(avg_train_values))\n",
    "    print(\"Training response range: \"+str(max_training_response_value-min_training_response_value))\n",
    "    logging.info(\"Training response range: \"+str(max_training_response_value-min_training_response_value))\n",
    "    print(\"R2 Training: \"+str(r2_training))\n",
    "    logging.info(\"R2 Training: \"+str(r2_training))\n",
    "    print(\"R2 Test (>0.7): \"+str(r2_test))\n",
    "    logging.info(\"R2 Test (>0.7): \"+str(r2_test))\n",
    "    print(\"LOO-Q2 (>0.7) : \"+str(float(LOO_Q2)))\n",
    "    logging.info(\"LOO-Q2 (>0.7) : \"+str(float(LOO_Q2)))\n",
    "    print(\"Q2-ext-F1 (>0.7): \"+str(q2_ext))\n",
    "    logging.info(\"Q2-ext-F1 (>0.7): \"+str(q2_ext))\n",
    "    print(\"SE Training: \"+str(se_training))\n",
    "    logging.info(\"SE Training: \"+str(se_training))\n",
    "    print(\"MAE Training: \"+str(mae_training))\n",
    "    logging.info(\"MAE Training: \"+str(mae_training))\n",
    "    print(\"MAE Test: \"+str(mae_test))\n",
    "    logging.info(\"MAE Test: \"+str(mae_test))\n",
    "    tropsha_pass = tropsha(test_value_inverse_transform_predicted,test_values_inverse_transform_actual)\n",
    "    if tropsha_pass:\n",
    "        print(\"Tropsha's Criteria: \"+\"Pass\")\n",
    "        logging.info(\"Tropsha's Criteria: \"+\"Pass\")\n",
    "    else:\n",
    "        print(\"Tropsha's Criteria: \"+\"Fail\")\n",
    "        logging.info(\"Tropsha's Criteria: \"+\"Fail\")\n",
    "    \n",
    "    # Save test outputs to evaluate MAE criteria. XternalValidationPlus can be used to evaluate MAE criteria.(https://sites.google.com/site/dtclabxvplus/)\n",
    "    data = [test_values_inverse_transform_actual]\n",
    "    data.append(test_value_inverse_transform_predicted)\n",
    "    data = np.transpose(data)\n",
    "    data_df = pd.DataFrame(data)\n",
    "    # For each fold of data Yobs and Ypred for test compounds is written to the following csv file. These files can be used to provide input to \n",
    "    # XternalValidationPlus (https://sites.google.com/site/dtclabxvplus/) to evaluate external validation parameters.\n",
    "    data_df.to_csv(str(out_folder_name)+\"/data_to_evaluate_xternal_validation_parameters/data_df\"+str(i)+\".csv\",header=['Yobs','Ypred'],index=False)\n",
    "\n",
    "    if (tropsha_pass and r2_test>0.7 and float(LOO_Q2)>0.7 and q2_ext>0.7):\n",
    "        \n",
    "        \n",
    "        print(\"Net: \"+\"Pass\")\n",
    "    else:\n",
    "        print(\"Net: \"+\"Fail\")\n",
    "        print(\"########\")\n",
    "        logging.info(\"########\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}