Machine_learning_models.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Author: Shayantan Banerjee\n",
    "# import warnings filter\n",
    "from warnings import simplefilter\n",
    "# ignore all future warnings\n",
    "simplefilter(action='ignore', category=FutureWarning)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from sklearn.preprocessing import Imputer\n",
    "import numpy as np\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import accuracy_score\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.feature_selection import SelectFromModel\n",
    "from sklearn.ensemble import ExtraTreesClassifier\n",
    "from sklearn.ensemble import GradientBoostingClassifier\n",
    "from sklearn.ensemble import AdaBoostClassifier\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.naive_bayes import MultinomialNB\n",
    "from sklearn.naive_bayes import GaussianNB\n",
    "from sklearn.svm import LinearSVC, SVC\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn.model_selection import cross_val_score\n",
    "from sklearn.model_selection import cross_validate\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.ensemble import RandomForestClassifier, VotingClassifier\n",
    "from sklearn import preprocessing\n",
    "from sklearn.preprocessing import MinMaxScaler\n",
    "import pickle\n",
    "from sklearn.metrics import confusion_matrix\n",
    "from sklearn.utils import resample\n",
    "from imblearn.over_sampling import SMOTE\n",
    "from sklearn.metrics import classification_report\n",
    "from sklearn import metrics\n",
    "from sklearn.metrics import roc_curve, auc\n",
    "from sklearn.metrics import roc_auc_score\n",
    "from sklearn.model_selection import cross_val_predict\n",
    "from sklearn.model_selection import StratifiedKFold\n",
    "import numpy as np\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn import datasets\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.feature_selection import SelectFromModel\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.feature_selection import RFECV\n",
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis\n",
    "from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis\n",
    "import scipy.stats as ss\n",
    "import itertools\n",
    "from sklearn.metrics import fbeta_score, make_scorer\n",
    "from sklearn.manifold import TSNE\n",
    "from sklearn.multiclass import OneVsRestClassifier\n",
    "from sklearn.model_selection import LeaveOneOut\n",
    "import mifs\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.model_selection import KFold\n",
    "import scipy\n",
    "import numpy as np, scipy.stats as st\n",
    "from sklearn.model_selection import RepeatedStratifiedKFold\n",
    "from sklearn.neural_network import MLPClassifier\n",
    "from scipy import interp\n",
    "from collections import Counter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "#importing list of DEGs from the gene expresssion analysis having pvalue <0.05 and |logfc| > 1\n",
    "deg1=pd.read_table(\"DEG1.0txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "#importing list of genes selected using L1 logistic based feature selection \n",
    "l1_reg=['AATBC', 'ACP6', 'ADGRB3', 'AMPD1', 'ANK3', 'ANKRD36C', 'AOC2', 'AREG',\n",
    "       'ARHGAP11B|LOC100288637', 'ARHGEF40', 'C12orf56', 'C19orf33', 'C1orf87',\n",
    "       'C6orf226', 'CAPZA3', 'CCL18', 'COL4A3', 'CPLANE1', 'CTC1', 'CTH',\n",
    "       'CTSL', 'CX3CR1', 'DCLK3', 'DCTN1-AS1', 'DDIT4', 'DDIT4L', 'DSCAML1',\n",
    "       'EPPK1', 'FAM222A-AS1', 'FAM66C', 'FAM92A', 'GALNT6', 'GPNMB', 'GPR171',\n",
    "       'GSTO1', 'GTPBP10', 'HIST1H3E', 'JUN', 'KIAA0586', 'KLHL23', 'KMO',\n",
    "       'LILRA4', 'LINC00518', 'LINC00893', 'LINC01137', 'LINC02433',\n",
    "       'LINC02544', 'LINC02579', 'LMO7DN', 'LOC100507516', 'LOC101927391',\n",
    "       'LOC101928487', 'LOC105376037', 'LOC283038', 'LOC284865', 'LOC401442',\n",
    "       'LRP2BP', 'MAGEB3', 'MAP3K7CL', 'MMP13', 'NCKAP5L', 'NECAB1', 'NLRP1',\n",
    "       'NXN|RHBDL2', 'OR7E24', 'PAX8-AS1', 'PCDHB15', 'PLCB1', 'PLTP',\n",
    "       'PLXNA4', 'PRAC1', 'PYGM', 'REG1B', 'RNF217-AS1', 'RTN1', 'SDC4',\n",
    "       'SEPT7P2', 'SLC25A43', 'SLC44A3', 'TCEAL9', 'TEAD3', 'TMEM176B',\n",
    "       'TNFAIP3', 'TTYH2', 'ZC2HC1B', 'ZNF484', 'ZNF674', 'ZXDA|ZXDB']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "#importing list of genes selected using MRMR based feature selection \n",
    "\n",
    "mrmr=['LRRC37A|LRRC37A3', 'PAM16', 'SDC4', 'PKHD1', 'KIF11', 'ULK2', 'BEND6',\n",
    "       'DAOA-AS1', 'HMHB1', 'MAS1L', 'TCN1', 'MAP3K7CL', 'SLC5A1', 'CAMK1D',\n",
    "       'WRAP53', 'SPRED2', 'SIGLEC16']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "#concatenating the three feature lists\n",
    "def Union(lst1, lst2): \n",
    "    final_list = lst1 + lst2 \n",
    "    return final_list\n",
    "feature_aggre=Union(Union(deg1['x'].tolist(),mrmr),l1_reg)\n",
    "mylist = list(dict.fromkeys(feature_aggre))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "#importing the gene expression values for 21724 genes\n",
    "data=pd.read_table(\"Gene_expression_without_controls.txt\")\n",
    "lab=pd.read_csv('labels_without_controls.txt')\n",
    "y=lab['x'].astype('category').cat.codes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "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>A1BG</th>\n",
       "      <th>A1BG-AS1</th>\n",
       "      <th>A1CF</th>\n",
       "      <th>A2M</th>\n",
       "      <th>A2M-AS1</th>\n",
       "      <th>A2ML1</th>\n",
       "      <th>A2MP1</th>\n",
       "      <th>A4GALT</th>\n",
       "      <th>A4GNT</th>\n",
       "      <th>AA06</th>\n",
       "      <th>...</th>\n",
       "      <th>ZWINT</th>\n",
       "      <th>ZXDA</th>\n",
       "      <th>ZXDA|ZXDB</th>\n",
       "      <th>ZXDB</th>\n",
       "      <th>ZXDC</th>\n",
       "      <th>ZYG11A</th>\n",
       "      <th>ZYG11B</th>\n",
       "      <th>ZYX</th>\n",
       "      <th>ZZEF1</th>\n",
       "      <th>ZZZ3</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>5.4783</td>\n",
       "      <td>5.5562</td>\n",
       "      <td>4.38005</td>\n",
       "      <td>5.36375</td>\n",
       "      <td>4.6992</td>\n",
       "      <td>3.78330</td>\n",
       "      <td>6.3970</td>\n",
       "      <td>7.1905</td>\n",
       "      <td>4.3406</td>\n",
       "      <td>5.0975</td>\n",
       "      <td>...</td>\n",
       "      <td>6.6838</td>\n",
       "      <td>4.6482</td>\n",
       "      <td>4.25410</td>\n",
       "      <td>4.43355</td>\n",
       "      <td>6.232683</td>\n",
       "      <td>5.0827</td>\n",
       "      <td>7.113900</td>\n",
       "      <td>9.28435</td>\n",
       "      <td>5.915833</td>\n",
       "      <td>5.96300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>5.5218</td>\n",
       "      <td>5.6292</td>\n",
       "      <td>4.45480</td>\n",
       "      <td>4.79720</td>\n",
       "      <td>4.6839</td>\n",
       "      <td>4.28380</td>\n",
       "      <td>6.3345</td>\n",
       "      <td>6.7634</td>\n",
       "      <td>4.0578</td>\n",
       "      <td>5.1819</td>\n",
       "      <td>...</td>\n",
       "      <td>7.7519</td>\n",
       "      <td>4.8541</td>\n",
       "      <td>4.22935</td>\n",
       "      <td>3.85635</td>\n",
       "      <td>6.228167</td>\n",
       "      <td>5.0560</td>\n",
       "      <td>8.274033</td>\n",
       "      <td>10.72825</td>\n",
       "      <td>6.008033</td>\n",
       "      <td>6.53555</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>5.3961</td>\n",
       "      <td>5.9805</td>\n",
       "      <td>4.80410</td>\n",
       "      <td>5.57900</td>\n",
       "      <td>5.6614</td>\n",
       "      <td>4.23930</td>\n",
       "      <td>6.9344</td>\n",
       "      <td>7.6834</td>\n",
       "      <td>4.3675</td>\n",
       "      <td>5.7251</td>\n",
       "      <td>...</td>\n",
       "      <td>6.2407</td>\n",
       "      <td>4.5151</td>\n",
       "      <td>4.25240</td>\n",
       "      <td>3.74830</td>\n",
       "      <td>5.457817</td>\n",
       "      <td>5.2001</td>\n",
       "      <td>6.370767</td>\n",
       "      <td>8.29315</td>\n",
       "      <td>5.708933</td>\n",
       "      <td>5.34015</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>5.4898</td>\n",
       "      <td>5.7417</td>\n",
       "      <td>4.56785</td>\n",
       "      <td>4.77960</td>\n",
       "      <td>5.5592</td>\n",
       "      <td>4.22800</td>\n",
       "      <td>5.4526</td>\n",
       "      <td>7.1383</td>\n",
       "      <td>4.1215</td>\n",
       "      <td>5.2230</td>\n",
       "      <td>...</td>\n",
       "      <td>8.6021</td>\n",
       "      <td>4.7943</td>\n",
       "      <td>4.16995</td>\n",
       "      <td>4.09635</td>\n",
       "      <td>6.235833</td>\n",
       "      <td>5.3694</td>\n",
       "      <td>7.648533</td>\n",
       "      <td>10.35275</td>\n",
       "      <td>5.583700</td>\n",
       "      <td>5.98235</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5.2333</td>\n",
       "      <td>5.6149</td>\n",
       "      <td>4.28090</td>\n",
       "      <td>5.50210</td>\n",
       "      <td>4.3983</td>\n",
       "      <td>3.96675</td>\n",
       "      <td>5.8352</td>\n",
       "      <td>6.8692</td>\n",
       "      <td>4.3052</td>\n",
       "      <td>5.0453</td>\n",
       "      <td>...</td>\n",
       "      <td>9.4406</td>\n",
       "      <td>5.0753</td>\n",
       "      <td>3.99805</td>\n",
       "      <td>4.51725</td>\n",
       "      <td>5.820800</td>\n",
       "      <td>4.7759</td>\n",
       "      <td>6.609800</td>\n",
       "      <td>9.96975</td>\n",
       "      <td>5.824167</td>\n",
       "      <td>6.68905</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>5.1145</td>\n",
       "      <td>5.1255</td>\n",
       "      <td>3.98305</td>\n",
       "      <td>4.86580</td>\n",
       "      <td>5.8344</td>\n",
       "      <td>4.22245</td>\n",
       "      <td>6.0768</td>\n",
       "      <td>6.6950</td>\n",
       "      <td>4.4106</td>\n",
       "      <td>5.1940</td>\n",
       "      <td>...</td>\n",
       "      <td>6.1439</td>\n",
       "      <td>5.1411</td>\n",
       "      <td>4.25165</td>\n",
       "      <td>4.00945</td>\n",
       "      <td>6.254250</td>\n",
       "      <td>5.1328</td>\n",
       "      <td>8.582033</td>\n",
       "      <td>8.98740</td>\n",
       "      <td>5.667000</td>\n",
       "      <td>6.84560</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>5.8635</td>\n",
       "      <td>5.7155</td>\n",
       "      <td>4.57770</td>\n",
       "      <td>5.12280</td>\n",
       "      <td>6.3146</td>\n",
       "      <td>3.80760</td>\n",
       "      <td>6.2055</td>\n",
       "      <td>7.3145</td>\n",
       "      <td>4.3032</td>\n",
       "      <td>4.5546</td>\n",
       "      <td>...</td>\n",
       "      <td>9.0421</td>\n",
       "      <td>5.1097</td>\n",
       "      <td>4.18775</td>\n",
       "      <td>4.66530</td>\n",
       "      <td>5.783433</td>\n",
       "      <td>4.5143</td>\n",
       "      <td>6.285133</td>\n",
       "      <td>9.07805</td>\n",
       "      <td>5.845700</td>\n",
       "      <td>6.61700</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>6.2135</td>\n",
       "      <td>5.2366</td>\n",
       "      <td>5.31575</td>\n",
       "      <td>5.73190</td>\n",
       "      <td>5.9888</td>\n",
       "      <td>4.17960</td>\n",
       "      <td>5.3496</td>\n",
       "      <td>7.4584</td>\n",
       "      <td>4.7521</td>\n",
       "      <td>5.6987</td>\n",
       "      <td>...</td>\n",
       "      <td>7.9171</td>\n",
       "      <td>5.1695</td>\n",
       "      <td>4.41530</td>\n",
       "      <td>4.93095</td>\n",
       "      <td>5.836017</td>\n",
       "      <td>5.8393</td>\n",
       "      <td>6.715467</td>\n",
       "      <td>9.05510</td>\n",
       "      <td>5.594733</td>\n",
       "      <td>6.45725</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>6.1481</td>\n",
       "      <td>6.1432</td>\n",
       "      <td>4.78990</td>\n",
       "      <td>5.62275</td>\n",
       "      <td>6.4451</td>\n",
       "      <td>4.55265</td>\n",
       "      <td>6.2823</td>\n",
       "      <td>6.8510</td>\n",
       "      <td>4.9676</td>\n",
       "      <td>5.5955</td>\n",
       "      <td>...</td>\n",
       "      <td>7.4816</td>\n",
       "      <td>4.8315</td>\n",
       "      <td>4.51725</td>\n",
       "      <td>3.98485</td>\n",
       "      <td>5.683133</td>\n",
       "      <td>5.6468</td>\n",
       "      <td>6.497933</td>\n",
       "      <td>8.20265</td>\n",
       "      <td>5.764067</td>\n",
       "      <td>5.86195</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>5.1902</td>\n",
       "      <td>5.3532</td>\n",
       "      <td>4.01435</td>\n",
       "      <td>4.96340</td>\n",
       "      <td>4.4271</td>\n",
       "      <td>3.81760</td>\n",
       "      <td>5.5444</td>\n",
       "      <td>6.6491</td>\n",
       "      <td>4.2248</td>\n",
       "      <td>5.0732</td>\n",
       "      <td>...</td>\n",
       "      <td>5.8829</td>\n",
       "      <td>4.2650</td>\n",
       "      <td>3.96990</td>\n",
       "      <td>3.99575</td>\n",
       "      <td>6.567100</td>\n",
       "      <td>4.8117</td>\n",
       "      <td>8.475333</td>\n",
       "      <td>10.97400</td>\n",
       "      <td>6.172133</td>\n",
       "      <td>6.46855</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>10 rows × 21724 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     A1BG  A1BG-AS1     A1CF      A2M  A2M-AS1    A2ML1   A2MP1  A4GALT  \\\n",
       "0  5.4783    5.5562  4.38005  5.36375   4.6992  3.78330  6.3970  7.1905   \n",
       "1  5.5218    5.6292  4.45480  4.79720   4.6839  4.28380  6.3345  6.7634   \n",
       "2  5.3961    5.9805  4.80410  5.57900   5.6614  4.23930  6.9344  7.6834   \n",
       "3  5.4898    5.7417  4.56785  4.77960   5.5592  4.22800  5.4526  7.1383   \n",
       "4  5.2333    5.6149  4.28090  5.50210   4.3983  3.96675  5.8352  6.8692   \n",
       "5  5.1145    5.1255  3.98305  4.86580   5.8344  4.22245  6.0768  6.6950   \n",
       "6  5.8635    5.7155  4.57770  5.12280   6.3146  3.80760  6.2055  7.3145   \n",
       "7  6.2135    5.2366  5.31575  5.73190   5.9888  4.17960  5.3496  7.4584   \n",
       "8  6.1481    6.1432  4.78990  5.62275   6.4451  4.55265  6.2823  6.8510   \n",
       "9  5.1902    5.3532  4.01435  4.96340   4.4271  3.81760  5.5444  6.6491   \n",
       "\n",
       "    A4GNT    AA06   ...      ZWINT    ZXDA  ZXDA|ZXDB     ZXDB      ZXDC  \\\n",
       "0  4.3406  5.0975   ...     6.6838  4.6482    4.25410  4.43355  6.232683   \n",
       "1  4.0578  5.1819   ...     7.7519  4.8541    4.22935  3.85635  6.228167   \n",
       "2  4.3675  5.7251   ...     6.2407  4.5151    4.25240  3.74830  5.457817   \n",
       "3  4.1215  5.2230   ...     8.6021  4.7943    4.16995  4.09635  6.235833   \n",
       "4  4.3052  5.0453   ...     9.4406  5.0753    3.99805  4.51725  5.820800   \n",
       "5  4.4106  5.1940   ...     6.1439  5.1411    4.25165  4.00945  6.254250   \n",
       "6  4.3032  4.5546   ...     9.0421  5.1097    4.18775  4.66530  5.783433   \n",
       "7  4.7521  5.6987   ...     7.9171  5.1695    4.41530  4.93095  5.836017   \n",
       "8  4.9676  5.5955   ...     7.4816  4.8315    4.51725  3.98485  5.683133   \n",
       "9  4.2248  5.0732   ...     5.8829  4.2650    3.96990  3.99575  6.567100   \n",
       "\n",
       "   ZYG11A    ZYG11B       ZYX     ZZEF1     ZZZ3  \n",
       "0  5.0827  7.113900   9.28435  5.915833  5.96300  \n",
       "1  5.0560  8.274033  10.72825  6.008033  6.53555  \n",
       "2  5.2001  6.370767   8.29315  5.708933  5.34015  \n",
       "3  5.3694  7.648533  10.35275  5.583700  5.98235  \n",
       "4  4.7759  6.609800   9.96975  5.824167  6.68905  \n",
       "5  5.1328  8.582033   8.98740  5.667000  6.84560  \n",
       "6  4.5143  6.285133   9.07805  5.845700  6.61700  \n",
       "7  5.8393  6.715467   9.05510  5.594733  6.45725  \n",
       "8  5.6468  6.497933   8.20265  5.764067  5.86195  \n",
       "9  4.8117  8.475333  10.97400  6.172133  6.46855  \n",
       "\n",
       "[10 rows x 21724 columns]"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "subset=data[mylist]\n",
    "X=subset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(228, 120)"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "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>MMP8</th>\n",
       "      <th>OLFM4</th>\n",
       "      <th>CEACAM8</th>\n",
       "      <th>LCN2</th>\n",
       "      <th>RETN</th>\n",
       "      <th>IL1R2</th>\n",
       "      <th>LTF</th>\n",
       "      <th>TCN1</th>\n",
       "      <th>ELANE</th>\n",
       "      <th>CLEC5A</th>\n",
       "      <th>...</th>\n",
       "      <th>SLC44A3</th>\n",
       "      <th>TCEAL9</th>\n",
       "      <th>TEAD3</th>\n",
       "      <th>TMEM176B</th>\n",
       "      <th>TNFAIP3</th>\n",
       "      <th>TTYH2</th>\n",
       "      <th>ZC2HC1B</th>\n",
       "      <th>ZNF484</th>\n",
       "      <th>ZNF674</th>\n",
       "      <th>ZXDA|ZXDB</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>11.92515</td>\n",
       "      <td>9.8809</td>\n",
       "      <td>9.7849</td>\n",
       "      <td>12.1490</td>\n",
       "      <td>12.3596</td>\n",
       "      <td>13.68595</td>\n",
       "      <td>11.8759</td>\n",
       "      <td>9.9560</td>\n",
       "      <td>5.9019</td>\n",
       "      <td>8.9816</td>\n",
       "      <td>...</td>\n",
       "      <td>3.7684</td>\n",
       "      <td>5.0308</td>\n",
       "      <td>3.6185</td>\n",
       "      <td>6.62240</td>\n",
       "      <td>9.84210</td>\n",
       "      <td>5.15370</td>\n",
       "      <td>3.5722</td>\n",
       "      <td>4.0165</td>\n",
       "      <td>2.4708</td>\n",
       "      <td>4.25410</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>12.68455</td>\n",
       "      <td>11.6035</td>\n",
       "      <td>9.3088</td>\n",
       "      <td>12.5925</td>\n",
       "      <td>12.0684</td>\n",
       "      <td>13.58880</td>\n",
       "      <td>12.1626</td>\n",
       "      <td>10.5258</td>\n",
       "      <td>8.6122</td>\n",
       "      <td>8.9461</td>\n",
       "      <td>...</td>\n",
       "      <td>3.1388</td>\n",
       "      <td>3.7964</td>\n",
       "      <td>3.9276</td>\n",
       "      <td>6.59895</td>\n",
       "      <td>9.41105</td>\n",
       "      <td>5.75325</td>\n",
       "      <td>3.6877</td>\n",
       "      <td>4.0340</td>\n",
       "      <td>2.4721</td>\n",
       "      <td>4.22935</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>8.66905</td>\n",
       "      <td>11.2306</td>\n",
       "      <td>5.9323</td>\n",
       "      <td>9.9195</td>\n",
       "      <td>10.7973</td>\n",
       "      <td>13.30405</td>\n",
       "      <td>9.3591</td>\n",
       "      <td>7.5340</td>\n",
       "      <td>6.1853</td>\n",
       "      <td>9.1950</td>\n",
       "      <td>...</td>\n",
       "      <td>3.2830</td>\n",
       "      <td>3.4904</td>\n",
       "      <td>3.8792</td>\n",
       "      <td>8.57845</td>\n",
       "      <td>11.73205</td>\n",
       "      <td>5.46100</td>\n",
       "      <td>3.3063</td>\n",
       "      <td>4.0867</td>\n",
       "      <td>2.5347</td>\n",
       "      <td>4.25240</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3 rows × 120 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       MMP8    OLFM4  CEACAM8     LCN2     RETN     IL1R2      LTF     TCN1  \\\n",
       "0  11.92515   9.8809   9.7849  12.1490  12.3596  13.68595  11.8759   9.9560   \n",
       "1  12.68455  11.6035   9.3088  12.5925  12.0684  13.58880  12.1626  10.5258   \n",
       "2   8.66905  11.2306   5.9323   9.9195  10.7973  13.30405   9.3591   7.5340   \n",
       "\n",
       "    ELANE  CLEC5A    ...      SLC44A3  TCEAL9   TEAD3  TMEM176B   TNFAIP3  \\\n",
       "0  5.9019  8.9816    ...       3.7684  5.0308  3.6185   6.62240   9.84210   \n",
       "1  8.6122  8.9461    ...       3.1388  3.7964  3.9276   6.59895   9.41105   \n",
       "2  6.1853  9.1950    ...       3.2830  3.4904  3.8792   8.57845  11.73205   \n",
       "\n",
       "     TTYH2  ZC2HC1B  ZNF484  ZNF674  ZXDA|ZXDB  \n",
       "0  5.15370   3.5722  4.0165  2.4708    4.25410  \n",
       "1  5.75325   3.6877  4.0340  2.4721    4.22935  \n",
       "2  5.46100   3.3063  4.0867  2.5347    4.25240  \n",
       "\n",
       "[3 rows x 120 columns]"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "#scaling of the data for selecting the best features using GridSearchCV\n",
    "from sklearn import preprocessing\n",
    "\n",
    "x = X #returns a numpy array\n",
    "min_max_scaler = preprocessing.MinMaxScaler()\n",
    "x_scaled = min_max_scaler.fit_transform(x)\n",
    "X_scaled = pd.DataFrame(x_scaled)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best Penalty: l2\n",
      "Best C: 2.154434690031882\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/root/anaconda3/lib/python3.7/site-packages/sklearn/model_selection/_search.py:841: DeprecationWarning: The default of the `iid` parameter will change from True to False in version 0.22 and will be removed in 0.24. This will change numeric results when test-set sizes are unequal.\n",
      "  DeprecationWarning)\n"
     ]
    }
   ],
   "source": [
    "#feature selection using GridsearchCV\n",
    "#Logit\n",
    "logreg=LogisticRegression()\n",
    "penalty = ['l1', 'l2']\n",
    "# Create regularization hyperparameter space\n",
    "C = np.logspace(-3, 3, 10)\n",
    "# Create hyperparameter options\n",
    "hyperparameters = dict(C=C, penalty=penalty)\n",
    "clf_logit = GridSearchCV(logreg, hyperparameters, cv=5, verbose=0)\n",
    "best_model = clf_logit.fit(X_scaled,y)\n",
    "print('Best Penalty:', best_model.best_estimator_.get_params()['penalty'])\n",
    "print('Best C:', best_model.best_estimator_.get_params()['C'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best parameters set found on development set:\n",
      "\n",
      "{'C': 1, 'kernel': 'linear'}\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/root/anaconda3/lib/python3.7/site-packages/sklearn/model_selection/_search.py:841: DeprecationWarning: The default of the `iid` parameter will change from True to False in version 0.22 and will be removed in 0.24. This will change numeric results when test-set sizes are unequal.\n",
      "  DeprecationWarning)\n"
     ]
    }
   ],
   "source": [
    "#feature selection using GridsearchCV\n",
    "#SVM\n",
    "SVM=SVC(kernel=\"linear\")\n",
    "tuned_parameters = [{'kernel': ['linear'], 'C': [1, 10, 100, 1000]}]\n",
    "clf_svm = GridSearchCV(SVC(), tuned_parameters, cv=5)\n",
    "best_model_svm=clf_svm.fit(X_scaled,y)\n",
    "print(\"Best parameters set found on development set:\")\n",
    "print()\n",
    "print(clf_svm.best_params_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/root/anaconda3/lib/python3.7/site-packages/sklearn/model_selection/_search.py:841: DeprecationWarning: The default of the `iid` parameter will change from True to False in version 0.22 and will be removed in 0.24. This will change numeric results when test-set sizes are unequal.\n",
      "  DeprecationWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best parameters set found on development set:\n",
      "\n",
      "{'criterion': 'gini', 'max_depth': 8, 'max_features': 'sqrt', 'n_estimators': 300}\n"
     ]
    }
   ],
   "source": [
    "#feature selection using GridsearchCV\n",
    "#Random forest\n",
    "param_grid = { \n",
    "    'n_estimators': [100, 300, 500],\n",
    "    'max_features': ['auto', 'sqrt', 'log2'],\n",
    "    'max_depth' : [4,5,6,7,8],\n",
    "    'criterion' :['gini', 'entropy']\n",
    "}\n",
    "clf_rf = GridSearchCV(RandomForestClassifier(), param_grid, cv=5)\n",
    "best_model_rf=clf_rf.fit(X_scaled,y)\n",
    "print(\"Best parameters set found on development set:\")\n",
    "print()\n",
    "print(clf_rf.best_params_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best parameters set found on development set:\n",
      "\n",
      "{'max_depth': 1, 'min_samples_split': 190}\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/root/anaconda3/lib/python3.7/site-packages/sklearn/model_selection/_search.py:841: DeprecationWarning: The default of the `iid` parameter will change from True to False in version 0.22 and will be removed in 0.24. This will change numeric results when test-set sizes are unequal.\n",
      "  DeprecationWarning)\n"
     ]
    }
   ],
   "source": [
    "#feature selection using Decision trees\n",
    "#Decision trees\n",
    "parameters={'min_samples_split' : range(10,500,20),'max_depth': range(1,20,2)}\n",
    "clf_dt = GridSearchCV(DecisionTreeClassifier(), parameters, cv=5)\n",
    "best_model_dt=clf_dt.fit(X_scaled,y)\n",
    "print(\"Best parameters set found on development set:\")\n",
    "print()\n",
    "print(clf_dt.best_params_)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "classifiers = {\n",
    "                'Logit': LogisticRegression(penalty='l2',C=2.15),\n",
    "                'SVM': SVC(C= 1,kernel= 'linear',probability=True),\n",
    "                'RandomForest': RandomForestClassifier(criterion='gini',max_depth=8,max_features='sqrt',n_estimators=300),\n",
    "                'DecisionTRees': DecisionTreeClassifier(max_depth=1,min_samples_split=190),\n",
    "                'LinearDiscriminantAnalysis': LinearDiscriminantAnalysis(),\n",
    "                                    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Running the repeated CV for stablity analysis of our models. Five different classifiers have been used\n",
    "# import warnings filter\n",
    "from warnings import simplefilter\n",
    "# ignore all future warnings\n",
    "simplefilter(action='ignore', category=FutureWarning)\n",
    "for i in classifiers:\n",
    "    print(i)\n",
    "    model=classifiers[i]\n",
    "    cross_validated(i,model,X,y,10)\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Repeated CV 50 times and averaging the classification metrics for all 50 runs. Also printing the best set of features\n",
    "#selected across all folds \n",
    "def cross_validated(j,model,X,y,k):\n",
    "    i=0\n",
    "    rskf = RepeatedStratifiedKFold(n_splits=10, n_repeats=50,random_state=36851234)\n",
    "    precision_comp=[]\n",
    "    recall_comp=[]\n",
    "    f1_comp=[]\n",
    "    roc_comp=[]\n",
    "    precision_uncomp=[]\n",
    "    recall_uncomp=[]\n",
    "    f1_uncomp=[]\n",
    "    tprs=[]\n",
    "    base_fpr = np.linspace(0, 1, 101)\n",
    "    fold=[]\n",
    "    o=[]\n",
    "    opti=[]\n",
    "    for train_index,test_index in rskf.split(X,y):\n",
    "        X_train,X_test=X.iloc[train_index],X.iloc[test_index]\n",
    "        y_train,y_test=y.iloc[train_index],y.iloc[test_index]\n",
    "        X_train_reduced, f,o=feature_reduction_using_RFECV(model,X_train,y_train)\n",
    "        fold.append(f)\n",
    "        opti.append(o)\n",
    "        X_test_reduced=X_test[X_train_reduced.columns]\n",
    "        model.fit(X_train_reduced,y_train)\n",
    "        ktest_pred_y = model.predict(X_test_reduced)\n",
    "        y_score = model.predict_proba(X_test_reduced)\n",
    "        fpr, tpr, _ = roc_curve(y_test, y_score[:, 1])\n",
    "        ktest_rocauc = roc_auc_score(y_test, ktest_pred_y)\n",
    "        report=classification_report(y_test,ktest_pred_y,target_names=[\"Not complicated\",\"Complicated\"],output_dict=True)\n",
    "        precision_comp.append(report[\"Complicated\"][\"precision\"])\n",
    "        recall_comp.append(report[\"Complicated\"][\"recall\"])\n",
    "        f1_comp.append(report['Complicated'][\"f1-score\"])\n",
    "        roc_comp.append(ktest_rocauc)\n",
    "        precision_uncomp.append(report[\"Not complicated\"][\"precision\"])\n",
    "        recall_uncomp.append(report[\"Not complicated\"][\"recall\"])\n",
    "        f1_uncomp.append(report['Not complicated'][\"f1-score\"])\n",
    "        \n",
    "    print(\"Precision UnComp \",np.mean(precision_uncomp),\"\\n\",\"Conf interval\",mean_confidence_interval(precision_uncomp), \"\\n\",\"Recall Uncomp\",np.mean(recall_uncomp),\"\\n\",\"Conf interval\",mean_confidence_interval(recall_uncomp),\"\\n\", \"F1 uncomp\",np.mean(f1_uncomp), \"\\n\",\"Conf interval\",mean_confidence_interval(f1_uncomp),\"\\n\",)\n",
    "    print(\"Precision comp\",np.mean(precision_comp),\"\\n\",\"Conf interval\",mean_confidence_interval(precision_comp),\"\\n\", \"Recall Comp\",np.mean(recall_comp),\"\\n\",\"Conf interval\",mean_confidence_interval(recall_comp), \"\\n\",\"F1 comp\",np.mean(f1_comp),\"\\n\", \"Conf interval\",mean_confidence_interval(f1_comp),\"\\n\",\"ROC \",np.mean(roc_comp),\"\\n\",\"Conf interval\",mean_confidence_interval(roc_comp))\n",
    "    flat_list = [item for sublist in fold for item in sublist]\n",
    "    merged = list(itertools.chain(*flat_list))\n",
    "    x=Counter(merged)\n",
    "    print(x)\n",
    "    print(\"Median optimal number of features \",np.median(opti))\n",
    "\n",
    "  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "#printing the CI for each classification metric\n",
    "def mean_confidence_interval(data, confidence=0.95):\n",
    "    a = 1.0 * np.array(data)\n",
    "    n = len(a)\n",
    "    m, se = np.mean(a), scipy.stats.sem(a)\n",
    "    h = se * scipy.stats.t.ppf((1 + confidence) / 2., n-1)\n",
    "    return m, m-h, m+h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Recursive feature elmination for each fold of the cross validation\n",
    "def feature_reduction_using_RFECV(model,X,y):\n",
    "    feature_list=[]\n",
    "    np.random.seed(315)\n",
    "    \n",
    "    rfecv = RFECV(estimator=model, step=1, cv=StratifiedKFold(10,random_state=1),\n",
    "                  scoring='roc_auc')\n",
    "    rfecv.fit(X, y)\n",
    "    names=X.columns # the names of the genes selected \n",
    "    new_X=X[X.columns[rfecv.support_==True]]\n",
    "    #returns the feature reduced dataset, number of features selcted for this run, names of the features selcted for this run\n",
    "    return new_X,rfecv.n_features_,new_X.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}