tutorials/Tutorial_KE07_Regression_ML_Stats.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/kamerlinlab/KIF/blob/main/tutorials/Tutorial_KE07_Regression_ML_Stats.ipynb)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tutorial: Regression ML and Statistical Analysis on the R1 and R4 KE07s\n",
    "\n",
    "In this jupyter notebook we will use the model_building.py and stat_modelling.py modules to perform machine learning (ML) and statisical analysis (both regression) on a kemp eliminase enzyme (KE07). Our target variable is a continous value (hence regression) and is the value of the KE07's W50 Chi2 angle. \n",
    "\n",
    "This notebook will also cover all the pre- and post-processing steps requireds to prepare, analyse and visualise the results.\n",
    "\n",
    "The dataset used here is for the R1 and R4 KE07s and is the same data as what was used in the manuscript. \n",
    "\n",
    "<center><img src=\"https://raw.githubusercontent.com/kamerlinlab/KIF/main/tutorials/miscellaneous/ke07_banner.png\" style=\"width: 70%\" /></center>\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup\n",
    "\n",
    "Install and load the required modules and then download the dataset we'll be working on from google drive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: KIF in c:\\users\\rory crean\\desktop\\github\\key-interactions-finder (0.1.1)\n",
      "Requirement already satisfied: pandas in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from KIF) (1.3.4)\n",
      "Requirement already satisfied: numpy in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from KIF) (1.21.4)\n",
      "Requirement already satisfied: scikit-learn in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from KIF) (1.0.1)\n",
      "Requirement already satisfied: scipy in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from KIF) (1.7.2)\n",
      "Requirement already satisfied: xgboost in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from KIF) (1.7.2)\n",
      "Requirement already satisfied: catboost in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from KIF) (1.1.1)\n",
      "Requirement already satisfied: MDAnalysis in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from KIF) (2.4.2)\n",
      "Requirement already satisfied: gdown in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from KIF) (4.6.0)\n",
      "Requirement already satisfied: plotly in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from catboost->KIF) (5.4.0)\n",
      "Requirement already satisfied: matplotlib in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from catboost->KIF) (3.5.0)\n",
      "Requirement already satisfied: graphviz in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from catboost->KIF) (0.20.1)\n",
      "Requirement already satisfied: six in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from catboost->KIF) (1.16.0)\n",
      "Requirement already satisfied: python-dateutil>=2.7.3 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from pandas->KIF) (2.8.2)\n",
      "Requirement already satisfied: pytz>=2017.3 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from pandas->KIF) (2021.3)\n",
      "Requirement already satisfied: tqdm in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from gdown->KIF) (4.64.1)\n",
      "Requirement already satisfied: requests[socks] in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from gdown->KIF) (2.28.1)\n",
      "Requirement already satisfied: beautifulsoup4 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from gdown->KIF) (4.11.1)\n",
      "Requirement already satisfied: filelock in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from gdown->KIF) (3.8.2)\n",
      "Requirement already satisfied: threadpoolctl in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from MDAnalysis->KIF) (3.0.0)\n",
      "Requirement already satisfied: networkx>=2.0 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from MDAnalysis->KIF) (2.8.8)\n",
      "Requirement already satisfied: mmtf-python>=1.0.0 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from MDAnalysis->KIF) (1.1.3)\n",
      "Requirement already satisfied: GridDataFormats>=0.4.0 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from MDAnalysis->KIF) (1.0.1)\n",
      "Requirement already satisfied: gsd>=1.9.3 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from MDAnalysis->KIF) (2.7.0)\n",
      "Requirement already satisfied: biopython>=1.80 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from MDAnalysis->KIF) (1.80)\n",
      "Requirement already satisfied: fasteners in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from MDAnalysis->KIF) (0.18)\n",
      "Requirement already satisfied: joblib>=0.12 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from MDAnalysis->KIF) (1.1.0)\n",
      "Requirement already satisfied: packaging in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from MDAnalysis->KIF) (21.3)\n",
      "Requirement already satisfied: mrcfile in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from GridDataFormats>=0.4.0->MDAnalysis->KIF) (1.4.3)\n",
      "Requirement already satisfied: pyparsing>=2.2.1 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from matplotlib->catboost->KIF) (3.0.6)\n",
      "Requirement already satisfied: cycler>=0.10 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from matplotlib->catboost->KIF) (0.11.0)\n",
      "Requirement already satisfied: fonttools>=4.22.0 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from matplotlib->catboost->KIF) (4.28.1)\n",
      "Requirement already satisfied: pillow>=6.2.0 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from matplotlib->catboost->KIF) (8.4.0)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from matplotlib->catboost->KIF) (1.3.2)\n",
      "Requirement already satisfied: msgpack>=1.0.0 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from mmtf-python>=1.0.0->MDAnalysis->KIF) (1.0.4)\n",
      "Requirement already satisfied: colorama in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from tqdm->gdown->KIF) (0.4.5)\n",
      "Requirement already satisfied: soupsieve>1.2 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from beautifulsoup4->gdown->KIF) (2.3.2.post1)\n",
      "Requirement already satisfied: tenacity>=6.2.0 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from plotly->catboost->KIF) (8.0.1)\n",
      "Requirement already satisfied: charset-normalizer<3,>=2 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from requests[socks]->gdown->KIF) (2.1.1)\n",
      "Requirement already satisfied: urllib3<1.27,>=1.21.1 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from requests[socks]->gdown->KIF) (1.26.13)\n",
      "Requirement already satisfied: certifi>=2017.4.17 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from requests[socks]->gdown->KIF) (2022.6.15)\n",
      "Requirement already satisfied: idna<4,>=2.5 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from requests[socks]->gdown->KIF) (3.4)\n",
      "Requirement already satisfied: PySocks!=1.5.7,>=1.5.6 in c:\\anaconda3\\envs\\ml_conf_features_3.8\\lib\\site-packages (from requests[socks]->gdown->KIF) (1.7.1)\n",
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "%pip install KIF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "from key_interactions_finder import data_preperation\n",
    "from key_interactions_finder import stat_modelling\n",
    "from key_interactions_finder import model_building\n",
    "from key_interactions_finder import post_proccessing\n",
    "from key_interactions_finder import pymol_projections"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will first need to donwload the dataset from google drive. \n",
    "The tutorial data will be saved in the relative path defined by \"save_dir\" in the cell block below.\n",
    "\n",
    "You can change this as you see fit. If you want to use the current directory you can do:\n",
    "\n",
    "save_dir=\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Downloading...\n",
      "From: https://drive.google.com/uc?id=13SjTIbSjF4ai_-Fn1vJDvW6-nsPmCdbw\n",
      "To: c:\\Users\\Rory Crean\\Desktop\\Github\\key-interactions-finder\\tutorials\\tutorial_datasets\\tutorial_dataset.zip\n",
      "100%|██████████| 22.9M/22.9M [00:00<00:00, 33.2MB/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tutorial files were successfully downloaded and unzipped.\n"
     ]
    }
   ],
   "source": [
    "from key_interactions_finder.utils import download_prep_tutorial_dataset\n",
    "\n",
    "drive_url = r\"https://drive.google.com/file/d/13SjTIbSjF4ai_-Fn1vJDvW6-nsPmCdbw/view?usp=share_link\"\n",
    "save_dir = \"tutorial_datasets/\"\n",
    "\n",
    "download_prep_tutorial_dataset(drive_url=drive_url, save_dir=save_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Where all input data is stored. \n",
    "in_dir = save_dir + r\"KE07_Tutorial/Input_data/\"\n",
    "\n",
    "# The target variable's per frame values are stored here. \n",
    "target_file = in_dir + r\"R1_5d2w_1in10_Trp50_Chi2.dat\"\n",
    "\n",
    "# Path to the variable we will use to filter frames with (optional addition for this system). \n",
    "w50_chi1_file = in_dir + r\"R1_5d2w_1in10_Trp50_Chi1.dat\"\n",
    "\n",
    "# output folders\n",
    "stats_out_dir = save_dir + r\"KE07_Tutorial/KE07_stat_analysis\" \n",
    "ml_out_dir = save_dir + r\"KE07_Tutorial/KE07_ml_analysis\""
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preperation Step 1: Load the non-covalent interaction datasets\n",
    "\n",
    "The contact identification calculation was split into 4 blocks of different residues ranges. We will first need to load these blocks in and merge them. Luckly this is very easy with pandas. \n",
    "\n",
    "Note this data was generated using the script: \"identify_contacts.py\" which is provided with KIF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>1Ala 4Lys Hbond</th>\n",
       "      <th>1Ala 5Arg Other</th>\n",
       "      <th>1Ala 219Asp Other</th>\n",
       "      <th>1Ala 247Asn Hbond</th>\n",
       "      <th>2Leu 247Asn Hbond</th>\n",
       "      <th>3Ala 219Asp Other</th>\n",
       "      <th>3Ala 247Asn Other</th>\n",
       "      <th>3Ala 248Val Hydrophobic</th>\n",
       "      <th>3Ala 249Arg Hbond</th>\n",
       "      <th>4Lys 45Asp Saltbr</th>\n",
       "      <th>...</th>\n",
       "      <th>240Tyr 243Lys Other</th>\n",
       "      <th>240Tyr 244His Hbond</th>\n",
       "      <th>241Leu 244His Other</th>\n",
       "      <th>241Leu 245Gly Other</th>\n",
       "      <th>241Leu 246Val Hydrophobic</th>\n",
       "      <th>241Leu 248Val Hydrophobic</th>\n",
       "      <th>242Lys 245Gly Other</th>\n",
       "      <th>242Lys 246Val Other</th>\n",
       "      <th>242Lys 248Val Hbond</th>\n",
       "      <th>242Lys 250Leu Hbond</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.0009</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.5091</td>\n",
       "      <td>1.4173</td>\n",
       "      <td>0.1942</td>\n",
       "      <td>8.4905</td>\n",
       "      <td>0.2583</td>\n",
       "      <td>0.8390</td>\n",
       "      <td>0.9887</td>\n",
       "      <td>...</td>\n",
       "      <td>4.7718</td>\n",
       "      <td>13.4418</td>\n",
       "      <td>4.6074</td>\n",
       "      <td>1.8850</td>\n",
       "      <td>6.7067</td>\n",
       "      <td>0.6970</td>\n",
       "      <td>2.3419</td>\n",
       "      <td>1.1636</td>\n",
       "      <td>1.4529</td>\n",
       "      <td>5.2773</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.3099</td>\n",
       "      <td>0.6169</td>\n",
       "      <td>0.1086</td>\n",
       "      <td>8.7308</td>\n",
       "      <td>0.0250</td>\n",
       "      <td>0.1596</td>\n",
       "      <td>0.4652</td>\n",
       "      <td>...</td>\n",
       "      <td>5.9753</td>\n",
       "      <td>17.8924</td>\n",
       "      <td>5.3730</td>\n",
       "      <td>1.3674</td>\n",
       "      <td>7.9847</td>\n",
       "      <td>1.3530</td>\n",
       "      <td>2.4137</td>\n",
       "      <td>1.2101</td>\n",
       "      <td>0.8352</td>\n",
       "      <td>3.9203</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.1025</td>\n",
       "      <td>4.1349</td>\n",
       "      <td>0.1174</td>\n",
       "      <td>4.0820</td>\n",
       "      <td>0.0004</td>\n",
       "      <td>0.0207</td>\n",
       "      <td>0.6851</td>\n",
       "      <td>...</td>\n",
       "      <td>4.2700</td>\n",
       "      <td>11.3290</td>\n",
       "      <td>4.1186</td>\n",
       "      <td>2.4458</td>\n",
       "      <td>5.9295</td>\n",
       "      <td>2.4387</td>\n",
       "      <td>2.8741</td>\n",
       "      <td>0.1052</td>\n",
       "      <td>0.4285</td>\n",
       "      <td>3.5174</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3 rows × 839 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   1Ala 4Lys Hbond  1Ala 5Arg Other  1Ala 219Asp Other  1Ala 247Asn Hbond  \\\n",
       "0           0.0009              0.0                0.0             0.5091   \n",
       "1           0.0000              0.0                0.0             0.3099   \n",
       "2           0.0000              0.0                0.0             1.1025   \n",
       "\n",
       "   2Leu 247Asn Hbond  3Ala 219Asp Other  3Ala 247Asn Other  \\\n",
       "0             1.4173             0.1942             8.4905   \n",
       "1             0.6169             0.1086             8.7308   \n",
       "2             4.1349             0.1174             4.0820   \n",
       "\n",
       "   3Ala 248Val Hydrophobic  3Ala 249Arg Hbond  4Lys 45Asp Saltbr  ...  \\\n",
       "0                   0.2583             0.8390             0.9887  ...   \n",
       "1                   0.0250             0.1596             0.4652  ...   \n",
       "2                   0.0004             0.0207             0.6851  ...   \n",
       "\n",
       "   240Tyr 243Lys Other  240Tyr 244His Hbond  241Leu 244His Other  \\\n",
       "0               4.7718              13.4418               4.6074   \n",
       "1               5.9753              17.8924               5.3730   \n",
       "2               4.2700              11.3290               4.1186   \n",
       "\n",
       "   241Leu 245Gly Other  241Leu 246Val Hydrophobic  241Leu 248Val Hydrophobic  \\\n",
       "0               1.8850                     6.7067                     0.6970   \n",
       "1               1.3674                     7.9847                     1.3530   \n",
       "2               2.4458                     5.9295                     2.4387   \n",
       "\n",
       "   242Lys 245Gly Other  242Lys 246Val Other  242Lys 248Val Hbond  \\\n",
       "0               2.3419               1.1636               1.4529   \n",
       "1               2.4137               1.2101               0.8352   \n",
       "2               2.8741               0.1052               0.4285   \n",
       "\n",
       "   242Lys 250Leu Hbond  \n",
       "0               5.2773  \n",
       "1               3.9203  \n",
       "2               3.5174  \n",
       "\n",
       "[3 rows x 839 columns]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "input_files = [\"KE07_block1.csv\", \"KE07_block2.csv\", \"KE07_block3.csv\", \"KE07_block4.csv\"]\n",
    "dfs = []\n",
    "for file_name in input_files:\n",
    "    file_path = in_dir + file_name\n",
    "    df = pd.read_csv(file_path)\n",
    "    dfs.append(df)\n",
    "\n",
    "all_contacts_df = pd.concat(dfs, join='outer', axis=1)\n",
    "all_contacts_df.head(3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see we now have a dataframe with all the contacts found (839) identified and of length 10001, with matches with the number of frames in the trajectory. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10001, 839)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_contacts_df.shape"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will add one more row to our dataframe which corresponds to the W50 Chi1 angle. This will allow us to filter the dataframe to remove those frames that belong to  \"Conformation B\". The reasons for this is  described in full in the paper, but this is essentially to help focus the analysis on the differences between the \"A\" and \"C\" state. \n",
    "- Note the target variable is the W50 Chi**2** angle, and we will filter data on the Chi**1** angle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>W50Chi1</th>\n",
       "      <th>1Ala 4Lys Hbond</th>\n",
       "      <th>1Ala 5Arg Other</th>\n",
       "      <th>1Ala 219Asp Other</th>\n",
       "      <th>1Ala 247Asn Hbond</th>\n",
       "      <th>2Leu 247Asn Hbond</th>\n",
       "      <th>3Ala 219Asp Other</th>\n",
       "      <th>3Ala 247Asn Other</th>\n",
       "      <th>3Ala 248Val Hydrophobic</th>\n",
       "      <th>3Ala 249Arg Hbond</th>\n",
       "      <th>...</th>\n",
       "      <th>240Tyr 243Lys Other</th>\n",
       "      <th>240Tyr 244His Hbond</th>\n",
       "      <th>241Leu 244His Other</th>\n",
       "      <th>241Leu 245Gly Other</th>\n",
       "      <th>241Leu 246Val Hydrophobic</th>\n",
       "      <th>241Leu 248Val Hydrophobic</th>\n",
       "      <th>242Lys 245Gly Other</th>\n",
       "      <th>242Lys 246Val Other</th>\n",
       "      <th>242Lys 248Val Hbond</th>\n",
       "      <th>242Lys 250Leu Hbond</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>193.934</td>\n",
       "      <td>0.0009</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.5091</td>\n",
       "      <td>1.4173</td>\n",
       "      <td>0.1942</td>\n",
       "      <td>8.4905</td>\n",
       "      <td>0.2583</td>\n",
       "      <td>0.8390</td>\n",
       "      <td>...</td>\n",
       "      <td>4.7718</td>\n",
       "      <td>13.4418</td>\n",
       "      <td>4.6074</td>\n",
       "      <td>1.8850</td>\n",
       "      <td>6.7067</td>\n",
       "      <td>0.6970</td>\n",
       "      <td>2.3419</td>\n",
       "      <td>1.1636</td>\n",
       "      <td>1.4529</td>\n",
       "      <td>5.2773</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>188.050</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.3099</td>\n",
       "      <td>0.6169</td>\n",
       "      <td>0.1086</td>\n",
       "      <td>8.7308</td>\n",
       "      <td>0.0250</td>\n",
       "      <td>0.1596</td>\n",
       "      <td>...</td>\n",
       "      <td>5.9753</td>\n",
       "      <td>17.8924</td>\n",
       "      <td>5.3730</td>\n",
       "      <td>1.3674</td>\n",
       "      <td>7.9847</td>\n",
       "      <td>1.3530</td>\n",
       "      <td>2.4137</td>\n",
       "      <td>1.2101</td>\n",
       "      <td>0.8352</td>\n",
       "      <td>3.9203</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>184.593</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.1025</td>\n",
       "      <td>4.1349</td>\n",
       "      <td>0.1174</td>\n",
       "      <td>4.0820</td>\n",
       "      <td>0.0004</td>\n",
       "      <td>0.0207</td>\n",
       "      <td>...</td>\n",
       "      <td>4.2700</td>\n",
       "      <td>11.3290</td>\n",
       "      <td>4.1186</td>\n",
       "      <td>2.4458</td>\n",
       "      <td>5.9295</td>\n",
       "      <td>2.4387</td>\n",
       "      <td>2.8741</td>\n",
       "      <td>0.1052</td>\n",
       "      <td>0.4285</td>\n",
       "      <td>3.5174</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3 rows × 840 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   W50Chi1  1Ala 4Lys Hbond  1Ala 5Arg Other  1Ala 219Asp Other  \\\n",
       "0  193.934           0.0009              0.0                0.0   \n",
       "1  188.050           0.0000              0.0                0.0   \n",
       "2  184.593           0.0000              0.0                0.0   \n",
       "\n",
       "   1Ala 247Asn Hbond  2Leu 247Asn Hbond  3Ala 219Asp Other  3Ala 247Asn Other  \\\n",
       "0             0.5091             1.4173             0.1942             8.4905   \n",
       "1             0.3099             0.6169             0.1086             8.7308   \n",
       "2             1.1025             4.1349             0.1174             4.0820   \n",
       "\n",
       "   3Ala 248Val Hydrophobic  3Ala 249Arg Hbond  ...  240Tyr 243Lys Other  \\\n",
       "0                   0.2583             0.8390  ...               4.7718   \n",
       "1                   0.0250             0.1596  ...               5.9753   \n",
       "2                   0.0004             0.0207  ...               4.2700   \n",
       "\n",
       "   240Tyr 244His Hbond  241Leu 244His Other  241Leu 245Gly Other  \\\n",
       "0              13.4418               4.6074               1.8850   \n",
       "1              17.8924               5.3730               1.3674   \n",
       "2              11.3290               4.1186               2.4458   \n",
       "\n",
       "   241Leu 246Val Hydrophobic  241Leu 248Val Hydrophobic  242Lys 245Gly Other  \\\n",
       "0                     6.7067                     0.6970               2.3419   \n",
       "1                     7.9847                     1.3530               2.4137   \n",
       "2                     5.9295                     2.4387               2.8741   \n",
       "\n",
       "   242Lys 246Val Other  242Lys 248Val Hbond  242Lys 250Leu Hbond  \n",
       "0               1.1636               1.4529               5.2773  \n",
       "1               1.2101               0.8352               3.9203  \n",
       "2               0.1052               0.4285               3.5174  \n",
       "\n",
       "[3 rows x 840 columns]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "just_chi1_df = pd.read_csv(w50_chi1_file)\n",
    "just_chi1_df = just_chi1_df.set_axis([\"W50Chi1\"], axis=1)\n",
    "\n",
    "chi1_df = pd.concat([just_chi1_df, all_contacts_df], axis=1)\n",
    "chi1_df.head(3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preperation Step 2. Prepare the Dataset for calculations with the data_preperation.py module. \n",
    "\n",
    "In this step, we take our dataframe and merge our per frame target file to it.\n",
    "\n",
    "We can also optionally perform several forms of filtering on the non-covalent interactions identified to select what types of interactions we would like to study.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Your PyContact features and target variable have been succesufully merged.\n",
      "You can access this dataset through the class attribute: '.df_processed'.\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Target</th>\n",
       "      <th>W50Chi1</th>\n",
       "      <th>1Ala 4Lys Hbond</th>\n",
       "      <th>1Ala 5Arg Other</th>\n",
       "      <th>1Ala 219Asp Other</th>\n",
       "      <th>1Ala 247Asn Hbond</th>\n",
       "      <th>2Leu 247Asn Hbond</th>\n",
       "      <th>3Ala 219Asp Other</th>\n",
       "      <th>3Ala 247Asn Other</th>\n",
       "      <th>3Ala 248Val Hydrophobic</th>\n",
       "      <th>...</th>\n",
       "      <th>240Tyr 243Lys Other</th>\n",
       "      <th>240Tyr 244His Hbond</th>\n",
       "      <th>241Leu 244His Other</th>\n",
       "      <th>241Leu 245Gly Other</th>\n",
       "      <th>241Leu 246Val Hydrophobic</th>\n",
       "      <th>241Leu 248Val Hydrophobic</th>\n",
       "      <th>242Lys 245Gly Other</th>\n",
       "      <th>242Lys 246Val Other</th>\n",
       "      <th>242Lys 248Val Hbond</th>\n",
       "      <th>242Lys 250Leu Hbond</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>85.8483</td>\n",
       "      <td>193.934</td>\n",
       "      <td>0.0009</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.5091</td>\n",
       "      <td>1.4173</td>\n",
       "      <td>0.1942</td>\n",
       "      <td>8.4905</td>\n",
       "      <td>0.2583</td>\n",
       "      <td>...</td>\n",
       "      <td>4.7718</td>\n",
       "      <td>13.4418</td>\n",
       "      <td>4.6074</td>\n",
       "      <td>1.8850</td>\n",
       "      <td>6.7067</td>\n",
       "      <td>0.6970</td>\n",
       "      <td>2.3419</td>\n",
       "      <td>1.1636</td>\n",
       "      <td>1.4529</td>\n",
       "      <td>5.2773</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>93.7505</td>\n",
       "      <td>188.050</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.3099</td>\n",
       "      <td>0.6169</td>\n",
       "      <td>0.1086</td>\n",
       "      <td>8.7308</td>\n",
       "      <td>0.0250</td>\n",
       "      <td>...</td>\n",
       "      <td>5.9753</td>\n",
       "      <td>17.8924</td>\n",
       "      <td>5.3730</td>\n",
       "      <td>1.3674</td>\n",
       "      <td>7.9847</td>\n",
       "      <td>1.3530</td>\n",
       "      <td>2.4137</td>\n",
       "      <td>1.2101</td>\n",
       "      <td>0.8352</td>\n",
       "      <td>3.9203</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>102.0940</td>\n",
       "      <td>184.593</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.1025</td>\n",
       "      <td>4.1349</td>\n",
       "      <td>0.1174</td>\n",
       "      <td>4.0820</td>\n",
       "      <td>0.0004</td>\n",
       "      <td>...</td>\n",
       "      <td>4.2700</td>\n",
       "      <td>11.3290</td>\n",
       "      <td>4.1186</td>\n",
       "      <td>2.4458</td>\n",
       "      <td>5.9295</td>\n",
       "      <td>2.4387</td>\n",
       "      <td>2.8741</td>\n",
       "      <td>0.1052</td>\n",
       "      <td>0.4285</td>\n",
       "      <td>3.5174</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3 rows × 841 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     Target  W50Chi1  1Ala 4Lys Hbond  1Ala 5Arg Other  1Ala 219Asp Other  \\\n",
       "0   85.8483  193.934           0.0009              0.0                0.0   \n",
       "1   93.7505  188.050           0.0000              0.0                0.0   \n",
       "2  102.0940  184.593           0.0000              0.0                0.0   \n",
       "\n",
       "   1Ala 247Asn Hbond  2Leu 247Asn Hbond  3Ala 219Asp Other  3Ala 247Asn Other  \\\n",
       "0             0.5091             1.4173             0.1942             8.4905   \n",
       "1             0.3099             0.6169             0.1086             8.7308   \n",
       "2             1.1025             4.1349             0.1174             4.0820   \n",
       "\n",
       "   3Ala 248Val Hydrophobic  ...  240Tyr 243Lys Other  240Tyr 244His Hbond  \\\n",
       "0                   0.2583  ...               4.7718              13.4418   \n",
       "1                   0.0250  ...               5.9753              17.8924   \n",
       "2                   0.0004  ...               4.2700              11.3290   \n",
       "\n",
       "   241Leu 244His Other  241Leu 245Gly Other  241Leu 246Val Hydrophobic  \\\n",
       "0               4.6074               1.8850                     6.7067   \n",
       "1               5.3730               1.3674                     7.9847   \n",
       "2               4.1186               2.4458                     5.9295   \n",
       "\n",
       "   241Leu 248Val Hydrophobic  242Lys 245Gly Other  242Lys 246Val Other  \\\n",
       "0                     0.6970               2.3419               1.1636   \n",
       "1                     1.3530               2.4137               1.2101   \n",
       "2                     2.4387               2.8741               0.1052   \n",
       "\n",
       "   242Lys 248Val Hbond  242Lys 250Leu Hbond  \n",
       "0               1.4529               5.2773  \n",
       "1               0.8352               3.9203  \n",
       "2               0.4285               3.5174  \n",
       "\n",
       "[3 rows x 841 columns]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Generate supervised dataset instance as we have a target variable.  \n",
    "supervised_dataset = data_preperation.SupervisedFeatureData(\n",
    "    input_df=chi1_df,\n",
    "    target_file=target_file,\n",
    "    is_classification=False,\n",
    "    header_present=True \n",
    ")\n",
    "\n",
    "supervised_dataset.df_processed.head(3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Optional Feature Filtering\n",
    "\n",
    "In the above dataframe we have 1887 columns (so 989 features + 1 target). We can take all of these forward for the stastical analysis or we can perform some filtering in advance (the choice is yours). \n",
    "There are five built in filtering methods available to you to perform filtering:\n",
    "\n",
    "1. **filter_by_occupancy(min_occupancy)** - Remove features that have an %occupancy less than the provided cut-off. %Occupancy is the % of frames with a non 0 value, i.e. the interaction is present in that frame.\n",
    "\n",
    "2. **filter_by_interaction_type(interaction_types_included)** - Inteactions are defined as one of four possible types: (\"Hbond\", \"Saltbr\", \"Hydrophobic\", \"Other\"). You select the interactions you want to include.\n",
    "\n",
    "3. **filter_by_avg_strength(average_strength_cut_off)** - Filter by the per frame contact score/strength for each interaction. You can filter features by the average score. Values below the cut-off are removed. \n",
    "\n",
    "4. **filter_by_occupancy_by_class(min_occupancy)** - Special alternative to the the standard filter features by occupancy method. %occupancy is determined for each class (as opposed to whole dataset), meaning only observations from 1 class have to meet the cut-off to keep the feature. Only avaible to datasets with a categorical target variable (classification). \n",
    "\n",
    "\n",
    "Finally if at any point in time you want to reset any filtering you've already performed, you can use the following method: \n",
    "\n",
    "5. **reset_filtering()** \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of features before any filtering: 840\n",
      "Number of features after filtering by occupancy: 830\n",
      "Number of features after filtering by average interaction scores: 684\n"
     ]
    }
   ],
   "source": [
    "supervised_dataset.reset_filtering() \n",
    "print(f\"Number of features before any filtering: {len(supervised_dataset.df_processed.columns) - 1}\")\n",
    "\n",
    "# Features with a %occupancy of less than 25% are removed. \n",
    "supervised_dataset.filter_by_occupancy(min_occupancy=25)\n",
    "print(f\"Number of features after filtering by occupancy: {len(supervised_dataset.df_filtered.columns) - 1}\")\n",
    "\n",
    "# Features with an average interaction strength less than 0.5 will be removed. \n",
    "supervised_dataset.filter_by_avg_strength(\n",
    "    average_strength_cut_off=0.5,  \n",
    ")\n",
    "print(f\"Number of features after filtering by average interaction scores: {len(supervised_dataset.df_filtered.columns) - 1}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now if we look at the class attributes of our SupervisedFeatureData() instance (we called it: supervised_dataset) using the special \"\\_\\_dict__\" method we can see two dataframes we could use in the stastical analysis to follow. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['input_df', 'is_classification', 'target_file', 'header_present', 'df_processed', 'df_filtered'])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "supervised_dataset.__dict__.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "They are: \n",
    "- 'df_processed' - The unfiltered dataframe, 1700 features\n",
    "- 'df_filtered' - The filtered dataframe. Less than 1700 features. \n",
    "\n",
    "In the following sections we will use the filtered dataframe but either dataframe could be justified based on your goals. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As described above, we will filter frames to remove those that belong to conformation B, using the W50 chi1 angle to define this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rows before filtering by W50 Chi1: 10001\n",
      "Rows after filtering by W50 Chi1: 9281\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Target</th>\n",
       "      <th>1Ala 247Asn Hbond</th>\n",
       "      <th>2Leu 247Asn Hbond</th>\n",
       "      <th>3Ala 247Asn Other</th>\n",
       "      <th>3Ala 248Val Hydrophobic</th>\n",
       "      <th>3Ala 249Arg Hbond</th>\n",
       "      <th>4Lys 45Asp Saltbr</th>\n",
       "      <th>4Lys 214Phe Hbond</th>\n",
       "      <th>4Lys 218Ala Hbond</th>\n",
       "      <th>4Lys 219Asp Saltbr</th>\n",
       "      <th>...</th>\n",
       "      <th>240Tyr 243Lys Other</th>\n",
       "      <th>240Tyr 244His Hbond</th>\n",
       "      <th>241Leu 244His Other</th>\n",
       "      <th>241Leu 245Gly Other</th>\n",
       "      <th>241Leu 246Val Hydrophobic</th>\n",
       "      <th>241Leu 248Val Hydrophobic</th>\n",
       "      <th>242Lys 245Gly Other</th>\n",
       "      <th>242Lys 246Val Other</th>\n",
       "      <th>242Lys 248Val Hbond</th>\n",
       "      <th>242Lys 250Leu Hbond</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>85.8483</td>\n",
       "      <td>0.5091</td>\n",
       "      <td>1.4173</td>\n",
       "      <td>8.4905</td>\n",
       "      <td>0.2583</td>\n",
       "      <td>0.8390</td>\n",
       "      <td>0.9887</td>\n",
       "      <td>6.2291</td>\n",
       "      <td>3.6381</td>\n",
       "      <td>5.2280</td>\n",
       "      <td>...</td>\n",
       "      <td>4.7718</td>\n",
       "      <td>13.4418</td>\n",
       "      <td>4.6074</td>\n",
       "      <td>1.8850</td>\n",
       "      <td>6.7067</td>\n",
       "      <td>0.6970</td>\n",
       "      <td>2.3419</td>\n",
       "      <td>1.1636</td>\n",
       "      <td>1.4529</td>\n",
       "      <td>5.2773</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>93.7505</td>\n",
       "      <td>0.3099</td>\n",
       "      <td>0.6169</td>\n",
       "      <td>8.7308</td>\n",
       "      <td>0.0250</td>\n",
       "      <td>0.1596</td>\n",
       "      <td>0.4652</td>\n",
       "      <td>6.7166</td>\n",
       "      <td>3.5377</td>\n",
       "      <td>5.8784</td>\n",
       "      <td>...</td>\n",
       "      <td>5.9753</td>\n",
       "      <td>17.8924</td>\n",
       "      <td>5.3730</td>\n",
       "      <td>1.3674</td>\n",
       "      <td>7.9847</td>\n",
       "      <td>1.3530</td>\n",
       "      <td>2.4137</td>\n",
       "      <td>1.2101</td>\n",
       "      <td>0.8352</td>\n",
       "      <td>3.9203</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>102.0940</td>\n",
       "      <td>1.1025</td>\n",
       "      <td>4.1349</td>\n",
       "      <td>4.0820</td>\n",
       "      <td>0.0004</td>\n",
       "      <td>0.0207</td>\n",
       "      <td>0.6851</td>\n",
       "      <td>7.3353</td>\n",
       "      <td>3.8191</td>\n",
       "      <td>6.3520</td>\n",
       "      <td>...</td>\n",
       "      <td>4.2700</td>\n",
       "      <td>11.3290</td>\n",
       "      <td>4.1186</td>\n",
       "      <td>2.4458</td>\n",
       "      <td>5.9295</td>\n",
       "      <td>2.4387</td>\n",
       "      <td>2.8741</td>\n",
       "      <td>0.1052</td>\n",
       "      <td>0.4285</td>\n",
       "      <td>3.5174</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3 rows × 684 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     Target  1Ala 247Asn Hbond  2Leu 247Asn Hbond  3Ala 247Asn Other  \\\n",
       "0   85.8483             0.5091             1.4173             8.4905   \n",
       "1   93.7505             0.3099             0.6169             8.7308   \n",
       "2  102.0940             1.1025             4.1349             4.0820   \n",
       "\n",
       "   3Ala 248Val Hydrophobic  3Ala 249Arg Hbond  4Lys 45Asp Saltbr  \\\n",
       "0                   0.2583             0.8390             0.9887   \n",
       "1                   0.0250             0.1596             0.4652   \n",
       "2                   0.0004             0.0207             0.6851   \n",
       "\n",
       "   4Lys 214Phe Hbond  4Lys 218Ala Hbond  4Lys 219Asp Saltbr  ...  \\\n",
       "0             6.2291             3.6381              5.2280  ...   \n",
       "1             6.7166             3.5377              5.8784  ...   \n",
       "2             7.3353             3.8191              6.3520  ...   \n",
       "\n",
       "   240Tyr 243Lys Other  240Tyr 244His Hbond  241Leu 244His Other  \\\n",
       "0               4.7718              13.4418               4.6074   \n",
       "1               5.9753              17.8924               5.3730   \n",
       "2               4.2700              11.3290               4.1186   \n",
       "\n",
       "   241Leu 245Gly Other  241Leu 246Val Hydrophobic  241Leu 248Val Hydrophobic  \\\n",
       "0               1.8850                     6.7067                     0.6970   \n",
       "1               1.3674                     7.9847                     1.3530   \n",
       "2               2.4458                     5.9295                     2.4387   \n",
       "\n",
       "   242Lys 245Gly Other  242Lys 246Val Other  242Lys 248Val Hbond  \\\n",
       "0               2.3419               1.1636               1.4529   \n",
       "1               2.4137               1.2101               0.8352   \n",
       "2               2.8741               0.1052               0.4285   \n",
       "\n",
       "   242Lys 250Leu Hbond  \n",
       "0               5.2773  \n",
       "1               3.9203  \n",
       "2               3.5174  \n",
       "\n",
       "[3 rows x 684 columns]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Remove all frames with W50 Chi1: < 160 and > 240. \n",
    "chi1_chi2_df = supervised_dataset.df_filtered\n",
    "filtered_df = chi1_chi2_df[( (chi1_chi2_df.W50Chi1 > 160) & (chi1_chi2_df.W50Chi1 < 240) )]\n",
    "\n",
    "# Now remove \"W50Chi1\" as its jobs is done. \n",
    "df_ready = (filtered_df.drop(\"W50Chi1\", axis=1)).reset_index(drop=True) \n",
    "\n",
    "print(f\"Rows before filtering by W50 Chi1: {len(supervised_dataset.df_filtered)}\")\n",
    "print(f\"Rows after filtering by W50 Chi1: {len(df_ready)}\")\n",
    "df_ready.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be seen above, we have now filtered our dataframe as desired and removed the W50Chi1 column from the dataframe so it is now ready for the stats and ML analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis Time! \n",
    "\n",
    "Our dataset is now ready for either statistical analysis or ML or both! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 1.1. Perform Statistical Analysis with the stat_modelling.py module. \n",
    "\n",
    "Now we will perform the actual statistical modelling to compare the differences for each feature against the target variable. \n",
    "\n",
    "With this module, we can calculate two different metrics to evaluate how dependant each feature is on the target variables value. They are:\n",
    "\n",
    "1. The [mutual information](https://scikit-learn.org/stable/modules/generated/sklearn.feature_selection.mutual_info_regression.html) using the implementation available in Scikit-learn. The mutual information can capture any kind of dependancy/relationship between variables and score their dependancy.\n",
    "\n",
    "2. The [linear correlation](). I assume this does not need an introduction. \n",
    "\n",
    "In both cases, the higher the absolute value, the more dependant the feature is on the target variable. The mutual information has scores in the range 0 and 1, whilst the linear correlation scores are in the range -1 to 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now time for stat regression analysis\n",
    "stat_model = stat_modelling.RegressionStatModel(\n",
    "    dataset=df_ready, # The dataframe we just made. \n",
    "    out_dir=stats_out_dir,\n",
    "    interaction_types_included=[\"Hbond\", \"Saltbr\", \"Hydrophobic\", \"Other\"] \n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can determine the values for our two metrics, these wont take long to run (maybe 2 mins in total)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mutual information scores calculated.\n",
      "tutorial_datasets\\KE07_Tutorial\\KE07_stat_analysis\\Mutual_Information_Per_Feature_Scores.csv written to disk.\n",
      "You can also access these results via the class attribute: 'mutual_infos'.\n"
     ]
    }
   ],
   "source": [
    "stat_model.calc_mutual_info_to_target()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Linear correlations calculated.\n",
      "tutorial_datasets\\KE07_Tutorial\\KE07_stat_analysis\\Linear_Correlations_Per_Feature_Scores.csv written to disk.\n",
      "You can also access these results via the class attribute: 'linear_correlations'.\n"
     ]
    }
   ],
   "source": [
    "stat_model.calc_linear_correl_to_target()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# As printed above we can access the results from these calculations from the class instance's (we called it stat_model) attributes. \n",
    "mi_results = stat_model.mutual_infos\n",
    "lc_results = stat_model.linear_correlations\n",
    "# stat_model.__dict__.keys() # uncomment to see all attributes available. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 1.2. Work up the Statistical Analysis with the post_proccessing.py module. \n",
    "\n",
    "In this module we can convert the per feature scores to per residue scores, by summing (and then normalising) every per feature score that each residue is involved in. This can allow us to identify residues which seem to differ the most between each state. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First generate an instance of the class. \n",
    "post_proc = post_proccessing.StatRegressorPostProcessor(\n",
    "    stat_model=stat_model,\n",
    "    out_dir=stats_out_dir\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tutorial_datasets\\KE07_Tutorial\\KE07_stat_analysis\\Mutual_Information_Scores_Per_Residue.csv written to disk.\n",
      "tutorial_datasets\\KE07_Tutorial\\KE07_stat_analysis\\Linear_Correlation_Scores_Per_Residue.csv written to disk.\n"
     ]
    }
   ],
   "source": [
    "# Now we can run the get_per_res_scores() method, changing the stat_method accordingly.\n",
    "mi_per_res_scores = post_proc.get_per_res_scores(\n",
    "    stat_method=\"mutual_information\")\n",
    "\n",
    "lin_correl_per_res_scores = post_proc.get_per_res_scores(\n",
    "    stat_method=\"linear_correlation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The per feature and per residue scores are not just saved to disk but available as variables so you can analyse them within python using whatever graphing program you like. \n",
    "\n",
    "For inspiration feel free to take a look at the article or the tutorial \"Tutorial_PTP1B_Classification_ML_Stats.ipynb\". "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 1.3. Project the Results onto Protein Structures with the pymol_projections.py module. \n",
    " \n",
    "Naturally, we may want to visualise some of the results we have generated above onto a protein structure. We can take advantage of\n",
    "the functions provided in the pymol_projections.py module to do this. \n",
    "\n",
    "As the name suggests this will output [PyMOL](https://pymol.org/) compatible python scripts which can be run to represent the results\n",
    "at the: \n",
    "\n",
    "1. Per feature level. (Cylinders are drawn between each feature, with the cylinder radii marking how strong the relative difference is. \n",
    "2. Per residue level. The Carbon alpha of each residue will be depicted as a sphere, with the sphere radii depicting how strong the the relative difference is. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The file: tutorial_datasets\\KE07_Tutorial\\KE07_stat_analysis\\linear_correlation_Pymol_Per_Feature_Scores.py was written to disk.\n",
      "The file: tutorial_datasets\\KE07_Tutorial\\KE07_stat_analysis\\mutual_information_Pymol_Per_Feature_Scores.py was written to disk.\n"
     ]
    }
   ],
   "source": [
    "# Write PyMOL compatable scripts for the per feature results.\n",
    "# Simply swap between the two statistical methods as shown below. \n",
    "pymol_projections.project_pymol_top_features(\n",
    "    per_feature_scores=stat_model.linear_correlations,\n",
    "    model_name=\"linear_correlation\",\n",
    "    numb_features=100, # top features to project, set to any integer or \"all\" for all features. \n",
    "    out_dir=stats_out_dir\n",
    ")\n",
    "\n",
    "pymol_projections.project_pymol_top_features(\n",
    "    per_feature_scores=stat_model.mutual_infos,\n",
    "    model_name=\"mutual_information\",\n",
    "    numb_features=100, # top features to project, set to any integer or \"all\" for all features. \n",
    "    out_dir=stats_out_dir\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The file: tutorial_datasets\\KE07_Tutorial\\KE07_stat_analysis\\linear_correlation_Pymol_Per_Res_Scores.py was written to disk.\n",
      "The file: tutorial_datasets\\KE07_Tutorial\\KE07_stat_analysis\\mutual_information_Pymol_Per_Res_Scores.py was written to disk.\n"
     ]
    }
   ],
   "source": [
    "# Write PyMOL compatable scripts for the per residue results.\n",
    "# Simply swap between the two statistical methods as shown below. \n",
    "pymol_projections.project_pymol_per_res_scores(\n",
    "    per_res_scores=lin_correl_per_res_scores,\n",
    "    model_name=\"linear_correlation\",\n",
    "    out_dir=stats_out_dir\n",
    ")  \n",
    "\n",
    "pymol_projections.project_pymol_per_res_scores(\n",
    "    per_res_scores=mi_per_res_scores,\n",
    "    model_name=\"mutual_information\",\n",
    "    out_dir=stats_out_dir\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are complete with the stats module. Here is an example of the kind of figures you can make with the pymol projections generated:\n",
    "\n",
    "\n",
    "<center><img src=\"https://raw.githubusercontent.com/kamerlinlab/KIF/main/tutorials/miscellaneous/ke07_example_outputs.png\" style=\"width: 70%\" /></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 2.1 Perform Machine Learning (ML) with the model_building.py module. \n",
    "\n",
    "Now we will use ML to generate models trained on our target variable. \n",
    "\n",
    "With this module, we use the feature importance scores from each ML model to to evaluate how different/similar each feature is against the target variable.\n",
    "\n",
    "**In the paper we used three ensemble based regression models:**\n",
    "\n",
    "1. [Categorical Boosting](https://catboost.ai/) - (Refered to as: CatBoost)\n",
    "\n",
    "2. [Extreme Gradient Boosting](https://xgboost.readthedocs.io/en/stable/) - (Refered to as: XGBoost)\n",
    "\n",
    "3. [Random Forest](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html)  (Refered to as: Random_Forest)\n",
    "\n",
    "**In this tutorial, we will only use the CatBoost algorithim, to reduce the time required**\n",
    "\n",
    "In all cases, the higher the score, the more \"different\" the feature is when in the two different states.\n",
    "\n",
    "We can use the same dataframe that we used for the statistical analysis in Part 1.1 below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Below is a summary of the machine learning you have planned.\n",
      "You will use 5-fold cross validation and perform 3 repeats.\n",
      "You will use up to 684 features to build each model, with 85.0% of your data used for training the model, which is 7888 observations. \n",
      "15.0% of your data will be used for evaluating the best models produced by the 5-fold cross validation, which is 1393 observations.\n",
      "You have selected to build 1 machine learning model(s), with the following hyperparameters: \n",
      " \n",
      "A CatBoost model, with grid search parameters: \n",
      "{'iterations': [100]} \n",
      "\n",
      "If you're happy with the above, lets get model building!\n"
     ]
    }
   ],
   "source": [
    "# Instantiate the ClassificationModel class. \n",
    "# Clearly there are many parameters here, using your IDE you can hover over RegressionModel to see what each parameter does. \n",
    "ml_model = model_building.RegressionModel(\n",
    "    dataset=df_ready,\n",
    "    evaluation_split_ratio=0.15,\n",
    "    models_to_use=[\"CatBoost\"], # \"XGBoost\", \"Random_Forest\"] # You can add the other methods back in you want. \n",
    "    scaling_method=\"min_max\",\n",
    "    out_dir=ml_out_dir, \n",
    "    cross_validation_splits=5, \n",
    "    cross_validation_repeats=3,\n",
    "    search_approach=\"none\",\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can go ahead and build the models.\n",
    "\n",
    "We have one optional parameter in the command below which is to save the models generated. This can be useful if you ever want to come back and do the post-processing later.\n",
    "\n",
    "If you set this to true all the files required will be saved to a folder called \"temporary_files\" in your current working directory. \n",
    "\n",
    "With the current setup this calculation will not take long (maybe 3 mins to run on a standard laptop). However, you could perform a very exhaustive calculation using grid search CV (possible by changing the \"search_approach\" parameter in model_building.RegressionModel() ), in which case it might be useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model saved to disk at: temporary_files\\CatBoost_Model.pickle\n",
      "Model building complete, returning final results with train/test datasets to you.\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>model</th>\n",
       "      <th>best_params</th>\n",
       "      <th>best_score</th>\n",
       "      <th>best_standard_deviation</th>\n",
       "      <th>Time taken to build model (minutes)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>CatBoost</td>\n",
       "      <td>{'iterations': 100}</td>\n",
       "      <td>0.869809</td>\n",
       "      <td>0.007601</td>\n",
       "      <td>2.62</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      model          best_params  best_score  best_standard_deviation  \\\n",
       "0  CatBoost  {'iterations': 100}    0.869809                 0.007601   \n",
       "\n",
       "   Time taken to build model (minutes)  \n",
       "0                                 2.62  "
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ml_model.build_models(save_models=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the model now built, we can see how long it took to build and some metrics describing the regression error for the train and test sets. \n",
    "\n",
    "We can now evaluate the quality of the model on the validation dataset (also sometimes refered to as the hold-out set)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Model</th>\n",
       "      <th>Explained Variance</th>\n",
       "      <th>Mean Absolute Error</th>\n",
       "      <th>MSE</th>\n",
       "      <th>RMSE</th>\n",
       "      <th>Mean Squared Log Error</th>\n",
       "      <th>r squared</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>CatBoost</td>\n",
       "      <td>0.8758</td>\n",
       "      <td>17.6354</td>\n",
       "      <td>700.6846</td>\n",
       "      <td>26.4704</td>\n",
       "      <td>0.0401</td>\n",
       "      <td>0.8756</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      Model  Explained Variance  Mean Absolute Error       MSE     RMSE  \\\n",
       "0  CatBoost              0.8758              17.6354  700.6846  26.4704   \n",
       "\n",
       "   Mean Squared Log Error  r squared  \n",
       "0                  0.0401     0.8756  "
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reports = ml_model.evaluate_models()\n",
    "reports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The report produced is a dataframe with 6 regressions metrics to enable you to evaluate the quality of the model. If you had more built more than one model this dataframe would contain additional rows for each additional model built.\n",
    "\n",
    "The MSE and RMSE stands for the mean squarred error and the root mean squared error respectively.\n",
    "\n",
    "Personally, I think the Mean Absolute Error (MAE) and RMSE are very useful metrics as they have the same units as your target dataset, meaning it is quite easy to use them to think about how good your model is. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2.2. Work up the ML results with the post_proccessing.py module. \n",
    "\n",
    "In order to perform the analysis we will need to provide the models generated in Part 2.1 Shown below are the two possible ways to do this. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First we will make an instance of the SupervisedPostProcessor class.\n",
    "ml_post_proc = post_proccessing.SupervisedPostProcessor(\n",
    "    out_dir=ml_out_dir,\n",
    ")\n",
    "\n",
    "# Option 1 - Load models from the instance of the SupervisedModel class. \n",
    "ml_post_proc.load_models_from_instance(supervised_model=ml_model)\n",
    "\n",
    "# Option 2 - Load models from disk. \n",
    "# (Use If you've already run the model building, shut down the kernel and now want to post-process).\n",
    "#ml_post_proc.load_models_from_disk(models_to_use=[\"XGBoost\", \"CatBoost\", \"Random_Forest\"]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tutorial_datasets\\KE07_Tutorial\\KE07_ml_analysis\\CatBoost_Feature_Scores.csv written to disk.\n",
      "All per feature scores have now been saved to disk.\n"
     ]
    }
   ],
   "source": [
    "# After preparing the class we can now determine the feature scores for each model made.\n",
    "ml_post_proc.get_per_feature_scores()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tutorial_datasets\\KE07_Tutorial\\KE07_ml_analysis\\CatBoost_Per_Residue_Scores.csv written to disk.\n",
      "All per residue scores have now been saved to disk.\n"
     ]
    }
   ],
   "source": [
    "# We can also project these per feature scores onto the per-residue level. \n",
    "ml_post_proc.get_per_res_scores()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['out_dir', 'feat_names', 'best_models', 'all_per_feature_scores', 'all_per_residue_scores'])\n"
     ]
    }
   ],
   "source": [
    "print(ml_post_proc.__dict__.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we take a look at the class attributes we can see the per feature and per residue scores were not just saved to disk, but are also stored inside the class.\n",
    "\n",
    "This means you could easily analyse them within Python if you want. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_per_res_scores = ml_post_proc.all_per_residue_scores\n",
    "all_feature_scores = ml_post_proc.all_per_feature_scores "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For ideas on how to work up these results please see the manuscript or our other tutorial: \"Tutorial_PTP1B_Classification_ML_Stats.ipynb\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 2.3. Project the Results onto Protein Structures with the pymol_projections.py module. \n",
    " \n",
    "This section is essentially identical to 1.3, only that now we will output the ml results instead of the stats results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The file: tutorial_datasets\\KE07_Tutorial\\KE07_ml_analysis\\CatBoost_Pymol_Per_Res_Scores.py was written to disk.\n",
      "The file: tutorial_datasets\\KE07_Tutorial\\KE07_ml_analysis\\CatBoost_Pymol_Per_Feature_Scores.py was written to disk.\n"
     ]
    }
   ],
   "source": [
    "# Here you do not need to specify what model you would like to output results for, all will be outputted simultaneously.\n",
    "pymol_projections.project_multiple_per_res_scores(\n",
    "    all_per_res_scores=ml_post_proc.all_per_residue_scores,\n",
    "    out_dir=ml_out_dir\n",
    ")\n",
    "\n",
    "pymol_projections.project_multiple_per_feature_scores(\n",
    "    all_per_feature_scores=ml_post_proc.all_per_feature_scores,\n",
    "    numb_features=100, # top features to project, set to any integer or \"all\" for all features.  \n",
    "    out_dir=ml_out_dir\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Done!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ML_conf_features_3.8",
   "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.8.12"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "d5d41bf74c17b5a6e59a3d560c52f1b05d1a02d51345af3a6d64602d4eaded4e"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}