ConvMixer-SKCM.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb8e43f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from tensorflow.keras import layers\n",
    "from tensorflow import keras\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import tensorflow_addons as tfa\n",
    "import tensorflow as tf\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4a1dc81",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select the time-to-event time\n",
    "def selected_time(select_time, time, event):\n",
    "    if (time > select_time): return 0\n",
    "    else:\n",
    "        if (event == 1.0): return 1\n",
    "        else: return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ca4f56b",
   "metadata": {},
   "outputs": [],
   "source": [
    "survival_metric_geo='Pat_Died'\n",
    "survival_time_geo='Pat_Overall_Survival_Months'\n",
    "time_geo=36"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3afc3b4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_clinical=pd.read_hdf('D:/MyData/Pathway_Attention/data/SKCM/train_set.h5',key='clinical')\n",
    "test_clinical=pd.read_hdf('data/SKCM/test_set.h5',key='clinical')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0de3ec5",
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_clinical.drop_duplicates()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b7a38bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(test_clinical.drop_duplicates())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33565510",
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_categorical = geo_clinical[[survival_metric_geo, survival_time_geo]].apply(lambda row: selected_time(time_geo, row[survival_time_geo], row[survival_metric_geo]), axis=1)\n",
    "geo_categorical.dropna(inplace=True)\n",
    "geo_categorical.value_counts(normalize=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99277299",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_cate= test_clinical[[survival_metric_geo, survival_time_geo]].apply(lambda row: selected_time(time_geo, row[survival_time_geo], row[survival_metric_geo]), axis=1)\n",
    "test_cate.dropna(inplace=True)\n",
    "test_cate.value_counts(normalize=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b5e9500",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df=pd.DataFrame(geo_categorical)\n",
    "df['Sam_Name']=geo_categorical.index\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9596cbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "dft=pd.DataFrame(test_cate)\n",
    "dft['Sam_Name']=test_cate.index\n",
    "dft"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "391a88ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(dft)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "581deda8",
   "metadata": {},
   "outputs": [],
   "source": [
    "df=df.drop_duplicates(subset=['Sam_Name'],keep='first')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efe700fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "278dcead",
   "metadata": {},
   "outputs": [],
   "source": [
    "dft=dft.drop_duplicates(subset=['Sam_Name'],keep='first')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea67a565",
   "metadata": {},
   "outputs": [],
   "source": [
    "dft"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfaed54e",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinical=geo_clinical.loc[[s for s in df.index]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f81169a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinical=clinical.drop_duplicates()\n",
    "clinical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e1001c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_clinical=test_clinical.loc[[s for s in dft.index]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73879d66",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "test_clinical=test_clinical.drop_duplicates()\n",
    "test_clinical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e2671be",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "#load the image data\n",
    "dir_name = 'D:/MyData/Pathway_Attention/data/IMG/'\n",
    "\n",
    "X = np.array([np.load(dir_name+s+\".npy\") for s in df[0].index]) \n",
    "X_test=np.array([np.load(dir_name+s+\".npy\") for s in dft[0].index]) \n",
    "'''\n",
    "X=[]\n",
    "for s in df[0].index:\n",
    "    array0=np.load(dir_name+s+\".npy\")\n",
    "    array1=np.stack([array0, array0, array0], axis=-1)\n",
    "    X.append(array1.tolist())\n",
    "X=np.array(X)\n",
    "\n",
    "X_test=[]\n",
    "for s in dft[0].index:\n",
    "    array0=np.load(dir_name+s+\".npy\")\n",
    "    array1=np.stack([array0, array0, array0], axis=-1)\n",
    "    X_test.append(array1.tolist())\n",
    "X_test=np.array(X_test)\n",
    "\n",
    "print(X.shape)\n",
    "print(X_test.shape)\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eab323eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(X.shape)\n",
    "print(X_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97ff90c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "from sklearn.preprocessing import LabelEncoder\n",
    "y = LabelEncoder().fit_transform(df[0])\n",
    "y=y.reshape(y.shape[0],1)\n",
    "y_test=LabelEncoder().fit_transform(dft[0])\n",
    "y_test=y_test.reshape(y_test.shape[0],1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83f1652a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "#Split into training set and test set\n",
    "from sklearn.model_selection import train_test_split\n",
    "X_train,X_val,y_train,y_val=train_test_split(X,y,test_size=0.2,random_state=1)\n",
    "print(X_train.shape)\n",
    "print(X_val.shape)\n",
    "print(y_train.shape)\n",
    "print(y_val.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "773a2c7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "X_random_train=X_train.reshape(X_train.shape[0],49*49)\n",
    "X_i_train=[[i] for i in range(len(y_train))]\n",
    "#X_smote.shape\n",
    "from imblearn.over_sampling import SMOTE as sm\n",
    "\n",
    "\n",
    "X_ri_train, y_r_train = sm(random_state=321).fit_resample(X_i_train,y_train)\n",
    "X_r_train=X_random_train[[i[0] for i in X_ri_train]]\n",
    "\n",
    "\n",
    "X_r_train=X_r_train.reshape(X_r_train.shape[0],49,49,1)\n",
    "\n",
    "y_r_train=y_r_train.reshape(y_r_train.shape[0],1)\n",
    "print(X_r_train.shape)\n",
    "print(y_r_train.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab57b72f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def activation_block(x):\n",
    "    x = layers.Activation(\"gelu\")(x)\n",
    "    return layers.BatchNormalization()(x)\n",
    "\n",
    "\n",
    "def conv_stem(x, filters: int, patch_size=7):\n",
    "    x = layers.Conv2D(filters, kernel_size=patch_size, strides=patch_size)(x)\n",
    "    return activation_block(x)\n",
    "\n",
    "\n",
    "def conv_mixer_block(x, filters: int, kernel_size: int):\n",
    "    # Depthwise convolution.\n",
    "    x0 = x\n",
    "    x = layers.DepthwiseConv2D(kernel_size=kernel_size, padding=\"same\")(x)\n",
    "    x = layers.Add()([activation_block(x), x0])  # Residual.\n",
    "\n",
    "    # Pointwise convolution.\n",
    "    x = layers.Conv2D(filters, kernel_size=1)(x)\n",
    "    x = layers.Activation(\"gelu\")(x)\n",
    "    x=layers.BatchNormalization()(x)\n",
    "    \n",
    "    return x\n",
    "\n",
    "def create_mask(mask_height, mask_width):  \n",
    "    mask = np.ones((mask_height, mask_width), dtype=np.float32)  \n",
    "    mask[-1, :] = 0  \n",
    "    if mask_height > 1 and mask_width > 1:  \n",
    "        mask[-2, -2:] = 0  \n",
    "    return tf.convert_to_tensor(mask)  #create mask to skip the added zero patches\n",
    "\n",
    "'''\n",
    "def get_conv_mixer_343_7(\n",
    "    image_size=49, filters=343, depth=7, kernel_size=3, patch_size=7, num_classes=2\n",
    "):\n",
    "    \"\"\"ConvMixer-256/8: https://openreview.net/pdf?id=TVHS5Y4dNvM.\n",
    "    The hyperparameter values are taken from the paper.\n",
    "    \"\"\"\n",
    "    inputs = keras.Input((image_size, image_size, 3))\n",
    "    #x = layers.Rescaling(scale=1.0 / 255)(inputs)\n",
    "   \n",
    "    # Extract patch embeddings.\n",
    "    x = conv_stem(inputs, filters, patch_size)\n",
    "\n",
    "    # ConvMixer blocks.\n",
    "    for _ in range(depth):\n",
    "        x = conv_mixer_block(x, filters, kernel_size)\n",
    "\n",
    "    # Classification block.\n",
    "    x = layers.GlobalAvgPool2D()(x)\n",
    "    outputs = layers.Dense(num_classes, activation=\"softmax\")(x)\n",
    "    \n",
    "    \n",
    "\n",
    "    return keras.Model(inputs, outputs)\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdac2428",
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_height = 7  \n",
    "mask_width = 7\n",
    "import keras_tuner as kt\n",
    "def buildConvMixer(hp):\n",
    "    inputs = keras.Input((49, 49, 1))\n",
    "    #the conv stem, looking for extracting yhe patch embeddings\n",
    "    x = layers.Conv2D(hp.Choice('filters1',[16,32,64,80,96,128,256]), kernel_size=7, strides=7)(inputs)\n",
    "    x = layers.Activation('gelu')(x)\n",
    "    x=layers.BatchNormalization()(x)\n",
    "    \n",
    "    x = conv_mixer_block(x, hp.Choice('filters1',[16,32,64,80,96,128,256]), kernel_size=hp.Choice('ks',[3,5,7]))\n",
    "    x = conv_mixer_block(x, hp.Choice('filters1',[16,32,64,80,96,128,256]), kernel_size=hp.Choice('ks',[3,5,7]))\n",
    "    #x = conv_mixer_block(x, hp.Choice('filters1',[16,32,64,128,256]), kernel_size=hp.Choice('ks',[3,5,7]))\n",
    "    #x = conv_mixer_block(x, hp.Choice('filters1',[32,48,64,80,96,128,256]), kernel_size=hp.Choice('ks',[3,5,7,9]))\n",
    "    #x = conv_mixer_block(x, hp.Choice('filters1',[64,80,96,112,128,160,256]), kernel_size=hp.Choice('ks',[3,5,7]))\n",
    "    \n",
    "    #mask the added zero patches\n",
    "    mask_tensor = create_mask(mask_height, mask_width)  \n",
    "    mask_broadcast = tf.expand_dims(tf.expand_dims(mask_tensor, axis=0), axis=-1)  # shape:(1, height, width, 1) \n",
    "    #mask_broadcast = tf.cast(mask_broadcast, dtype=tf.float16)\n",
    "    x_masked_broadcast = x * mask_broadcast  # automatically broadcast\n",
    "\n",
    "    # Classification block.\n",
    "    #x = layers.GlobalAvgPool2D()(x)\n",
    "    x = layers.Flatten()(x_masked_broadcast)\n",
    "    outputs = layers.Dense(1, activation=\"sigmoid\")(x)\n",
    "\n",
    "    model=keras.Model(inputs,outputs)\n",
    "     \n",
    "    optimizer = tfa.optimizers.AdamW(learning_rate=0.001, weight_decay=0.0001)\n",
    "\n",
    "    model.compile(\n",
    "        optimizer=optimizer,\n",
    "        loss=keras.losses.BinaryCrossentropy(\n",
    "            from_logits=False, label_smoothing=0.1\n",
    "        ),\n",
    "        metrics=[\n",
    "#            keras.metrics.BinaryAccuracy(name=\"accuracy\"),\n",
    "#            keras.metrics.BinaryAccuracy(\n",
    "#                name=\"binary_accuracy\", dtype=None, threshold=0.5\n",
    "#            ),\n",
    "            keras.metrics.AUC(name='auc',from_logits=False),\n",
    "#            keras.metrics.Precision(name='precision'),\n",
    "#            keras.metrics.Recall(name='recall'),  \n",
    "#            tfa.metrics.F1Score(num_classes,threshold=0.5,name='f1_score')\n",
    "        ],\n",
    "    )\n",
    "    \n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01d31a9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "tuner1 = kt.BayesianOptimization(\n",
    "        buildConvMixer,\n",
    "        objective=kt.Objective(\"val_auc\", direction=\"max\"),\n",
    "        max_trials=20,\n",
    "        overwrite=True)\n",
    "stop_early = tf.keras.callbacks.EarlyStopping(monitor='val_auc',mode='max', patience=10)\n",
    "tuner1.search([X_r_train],\n",
    "              y_r_train,\n",
    "              epochs=50,\n",
    "              #batch_size=hp.Choice(\"batch_size\", [16, 32]),\n",
    "              validation_data=([X_val], y_val),\n",
    "              callbacks=[stop_early])\n",
    "best_model = tuner1.get_best_models()[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "874f3078",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_hyperparameters = tuner1.get_best_hyperparameters(1)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e6a9794",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d78d5079",
   "metadata": {},
   "outputs": [],
   "source": [
    "tuner1.results_summary(num_trials=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4694d492",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sklearn.metrics\n",
    "y_pred=best_model.predict([X_test])\n",
    "auc=sklearn.metrics.roc_auc_score(y_test,y_pred)\n",
    "auc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b4c4dc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(best_model.predict([X_test]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "888522a0",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "# Code reference: https://bit.ly/3awIRbP.\n",
    "import matplotlib.cm as c\n",
    "\n",
    "def visualization_plot(weights, idx=1):\n",
    "    # First, apply min-max normalization to the\n",
    "    # given weights to avoid isotrophic scaling.\n",
    "    p_min, p_max = weights.min(), weights.max()\n",
    "    weights = (weights - p_min) / (p_max - p_min)\n",
    "\n",
    "    # Visualize all the filters.\n",
    "    num_filters = 16\n",
    "    plt.figure(figsize=(14, 14))\n",
    "\n",
    "    for i in range(num_filters):\n",
    "        current_weight = weights[:, :, :, i]\n",
    "        if current_weight.shape[-1] == 1:\n",
    "            current_weight = current_weight.squeeze()\n",
    "        ax = plt.subplot(4, 4, idx)\n",
    "        ax.set_xticks([])\n",
    "        ax.set_yticks([])\n",
    "        plt.imshow(current_weight,c.RdGy)\n",
    "        \n",
    "        idx += 1\n",
    "    plt.savefig('D:/Mydata/figs/feature_embeddings.svg', dpi=600, format='svg')\n",
    "\n",
    "\n",
    "\n",
    "# We first visualize the learned patch embeddings.\n",
    "patch_embeddings = best_model.layers[1].get_weights()[0]\n",
    "visualization_plot(patch_embeddings)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "478168b1",
   "metadata": {},
   "source": [
    "## Cross Validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ab1be5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "from sklearn.model_selection import train_test_split\n",
    "from imblearn.over_sampling import SMOTE as os\n",
    "from sklearn.model_selection import StratifiedKFold\n",
    "from sklearn.metrics import roc_curve\n",
    "import sklearn.metrics\n",
    "accuracy_list=[]\n",
    "auc_list=[]\n",
    "precision_list=[]\n",
    "recall_list=[]\n",
    "f1_score_list=[]\n",
    "thresh_list=[]\n",
    "\n",
    "\n",
    "accuracy_list_t=[]\n",
    "auc_list_t=[]\n",
    "precision_list_t=[]\n",
    "recall_list_t=[]\n",
    "f1_score_list_t=[]\n",
    "thresh_list_t=[]\n",
    "\n",
    "init_random_state=1\n",
    "for i in range(0,50):\n",
    "    X_train,X_val,geo_categorical_train,geo_categorical_test,y_train,y_val=train_test_split(X,df[0],y,test_size=0.2,random_state=init_random_state+i)\n",
    "    print('=============================================================================================================')\n",
    "    print('* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *')\n",
    "    print(' * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ')\n",
    "    print('* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *')\n",
    "    print(' * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ')\n",
    "    print('* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *')\n",
    "    print('=============================================================================================================')\n",
    "    print('current random state is: {}'.format(init_random_state+i))\n",
    "    print('This is the {}th fold'.format(i+1))\n",
    "   \n",
    "    X_random_train=X_train.reshape(X_train.shape[0],49*49)\n",
    "    \n",
    "    X_i_train=[[i] for i in range(len(y_train))]\n",
    "    X_ri_train, y_r_train = os(random_state=321+i).fit_resample(X_i_train,y_train)\n",
    "    X_r_train=X_random_train[[i[0] for i in X_ri_train]]\n",
    "    \n",
    "\n",
    "    X_r_train=X_r_train.reshape(X_r_train.shape[0],49,49,1)\n",
    "    \n",
    "    y_r_train=y_r_train.reshape(y_r_train.shape[0],1)\n",
    "\n",
    "    model=buildConvMixer(best_hyperparameters)\n",
    "    \n",
    "    lr_schedule = keras.optimizers.schedules.ExponentialDecay(\n",
    "    initial_learning_rate=4e-4,\n",
    "    decay_steps=16,\n",
    "    decay_rate=0.99)\n",
    "\n",
    "    \n",
    "    optimizer = tfa.optimizers.AdamW(learning_rate=lr_schedule, weight_decay=0.0001)\n",
    "    #optimizer=keras.optimizers.Adam(learning_rate=0.001)\n",
    "    model.compile(\n",
    "        optimizer=optimizer,\n",
    "        loss=keras.losses.BinaryCrossentropy(\n",
    "            from_logits=False, label_smoothing=0.1\n",
    "        ),\n",
    "        metrics=[\n",
    "#            keras.metrics.BinaryAccuracy(name=\"accuracy\"),\n",
    "#            keras.metrics.BinaryAccuracy(\n",
    "#                name=\"binary_accuracy\", dtype=None, threshold=0.5\n",
    "#            ),\n",
    "            keras.metrics.AUC(name='auc',from_logits=False,num_thresholds=500),\n",
    "#            keras.metrics.Precision(name='precision'),\n",
    "#            keras.metrics.Recall(name='recall'),  \n",
    "#            tfa.metrics.F1Score(num_classes,threshold=0.5,name='f1_score')\n",
    "        ],\n",
    "    )\n",
    "\n",
    "    checkpoint_filepath = \"./tmp/checkpoint\"\n",
    "    checkpoint_callback = keras.callbacks.ModelCheckpoint(\n",
    "        checkpoint_filepath,\n",
    "        monitor=\"val_auc\",\n",
    "        mode='max',\n",
    "        save_best_only=True,\n",
    "        save_weights_only=True,\n",
    "    )\n",
    "    \n",
    "    early_stop = keras.callbacks.EarlyStopping(monitor='val_auc',mode='max', patience=30)\n",
    "\n",
    "    history = model.fit(\n",
    "        x=([X_r_train]),\n",
    "        y=y_r_train,\n",
    "        batch_size=32,\n",
    "        epochs=150,\n",
    "        validation_data=([X_val], y_val),\n",
    "        callbacks=[checkpoint_callback,early_stop],\n",
    "    )\n",
    "\n",
    "    model.load_weights(checkpoint_filepath)\n",
    "#    _, accuracy, binary_accuracy,auc_roc = model.evaluate(X_test, y_test)\n",
    "    #_, accuracy,auc,precision,recall = model.evaluate([X_f_test,X_p_test,age_test,stage_test], y_test)\n",
    "    y_pred=model.predict([X_val])\n",
    "    auc=sklearn.metrics.roc_auc_score(y_val,y_pred)\n",
    "    geo_categorical_pred=pd.DataFrame(geo_categorical_test)\n",
    "    geo_categorical_pred['prediction']=y_pred\n",
    "    \n",
    "    fpr, tpr, thresholds = roc_curve(y_val, y_pred)\n",
    "    J = tpr - fpr\n",
    "    index_max = np.argmax(J)\n",
    "    best_thresh = thresholds[index_max]\n",
    "    \n",
    "    tp=geo_categorical_pred[(geo_categorical_pred[0]==1.0) & (geo_categorical_pred['prediction']>=best_thresh)].shape[0]\n",
    "    fp=geo_categorical_pred[(geo_categorical_pred[0]==1.0) & (geo_categorical_pred['prediction']<best_thresh)].shape[0]\n",
    "    tn=geo_categorical_pred[(geo_categorical_pred[0]==0.0) & (geo_categorical_pred['prediction']<=best_thresh)].shape[0]\n",
    "    fn=geo_categorical_pred[(geo_categorical_pred[0]==0.0) & (geo_categorical_pred['prediction']>best_thresh)].shape[0]\n",
    "    \n",
    "    accuracy=(tp+tn)/(tp+fp+tn+fn)\n",
    "    precision=tp/(tp+fp)\n",
    "    recall=tp/(tp+fn)\n",
    "    f1_score=2*(precision*recall)/(precision+recall) \n",
    "    \n",
    "#External Validation / Generalization Validation   \n",
    "    print('External Validation, waiting...')\n",
    "    y_pred_t=model.predict([X_test])\n",
    "    auc_t=sklearn.metrics.roc_auc_score(y_test,y_pred_t)\n",
    "    df_pred=pd.DataFrame(dft[0])\n",
    "    df_pred['prediction']=y_pred_t\n",
    "    \n",
    "    fpr_t, tpr_t, thresholds_t = roc_curve(y_test, y_pred_t)\n",
    "    J_t = tpr_t - fpr_t\n",
    "    index_max_t = np.argmax(J_t)\n",
    "    best_thresh_t = thresholds_t[index_max_t]\n",
    "    \n",
    "    tp_t=df_pred[(df_pred[0]==1.0) & (df_pred['prediction']>=best_thresh_t)].shape[0]\n",
    "    fp_t=df_pred[(df_pred[0]==1.0) & (df_pred['prediction']<best_thresh_t)].shape[0]\n",
    "    tn_t=df_pred[(df_pred[0]==0.0) & (df_pred['prediction']<=best_thresh_t)].shape[0]\n",
    "    fn_t=df_pred[(df_pred[0]==0.0) & (df_pred['prediction']>best_thresh_t)].shape[0]\n",
    "    \n",
    "    accuracy_t=(tp_t+tn_t)/(tp_t+fp_t+tn_t+fn_t)\n",
    "    precision_t=tp_t/(tp_t+fp_t)\n",
    "    recall_t=tp_t/(tp_t+fn_t)\n",
    "    f1_score_t=2*(precision_t*recall_t)/(precision_t+recall_t) \n",
    "    '''\n",
    "    y_pred_geo=model.predict([X_f_geo])\n",
    "    auc_geo=sklearn.metrics.roc_auc_score(y_geo,y_pred_geo)\n",
    "    geo_categorical_pred=pd.DataFrame(geo_categorical)\n",
    "    geo_categorical_pred['prediction']=y_pred_geo\n",
    "    \n",
    "    fpr_geo, tpr_geo, thresholds_geo = roc_curve(y_geo, y_pred_geo)\n",
    "    J_geo = tpr_geo - fpr_geo\n",
    "    index_max_geo = np.argmax(J_geo)\n",
    "    best_thresh_geo = thresholds_geo[index_max_geo]\n",
    "    \n",
    "    tp_geo=geo_categorical_pred[(geo_categorical_pred[0]==1.0) & (geo_categorical_pred['prediction']>=best_thresh_geo)].shape[0]\n",
    "    fp_geo=geo_categorical_pred[(geo_categorical_pred[0]==1.0) & (geo_categorical_pred['prediction']<best_thresh_geo)].shape[0]\n",
    "    tn_geo=geo_categorical_pred[(geo_categorical_pred[0]==0.0) & (geo_categorical_pred['prediction']<=best_thresh_geo)].shape[0]\n",
    "    fn_geo=geo_categorical_pred[(geo_categorical_pred[0]==0.0) & (geo_categorical_pred['prediction']>best_thresh_geo)].shape[0]\n",
    "    \n",
    "    accuracy_geo=(tp_geo+tn_geo)/(tp_geo+fp_geo+tn_geo+fn_geo)\n",
    "    precision_geo=tp_geo/(tp_geo+fp_geo)\n",
    "    recall_geo=tp_geo/(tp_geo+fn_geo)\n",
    "    f1_score_geo=2*(precision_geo*recall_geo)/(precision_geo+recall_geo) \n",
    "    '''\n",
    "   \n",
    "    #print(f\"Test accuracy: {round(accuracy * 100, 2)}%\")\n",
    "#    print(f\"Test binary accuracy: {round(binary_accuracy * 100, 2)}%\")\n",
    "    #print(f\"Test AUC: {round(auc, 4)}\")\n",
    "    #print(f\"Test precision: {round(precision * 100, 2)}%\")\n",
    "    #print(f\"Test recall: {round(recall * 100, 2)}%\")\n",
    "    #print(f\"Test f1_score: {round(f1_score*100,2)}%\")\n",
    "    accuracy_list.append(accuracy)\n",
    "    auc_list.append(auc)\n",
    "    precision_list.append(precision)\n",
    "    recall_list.append(recall)\n",
    "    f1_score_list.append(f1_score)\n",
    "    thresh_list.append(best_thresh)\n",
    "    \n",
    "    accuracy_list_t.append(accuracy_t)\n",
    "    auc_list_t.append(auc_t)\n",
    "    precision_list_t.append(precision_t)\n",
    "    recall_list_t.append(recall_t)\n",
    "    f1_score_list_t.append(f1_score_t)\n",
    "    thresh_list_t.append(best_thresh_t)\n",
    "    '''\n",
    "    accuracy_list_geo.append(accuracy_geo)\n",
    "    auc_list_geo.append(auc_geo)\n",
    "    precision_list_geo.append(precision_geo)\n",
    "    recall_list_geo.append(recall_geo)\n",
    "    f1_score_list_geo.append(f1_score_geo)\n",
    "    thresh_list_geo.append(best_thresh_geo)\n",
    "    '''\n",
    "print('mean accuracy={}±{}'.format(np.mean(accuracy_list),np.std(accuracy_list)))\n",
    "print('mean auc={}±{}'.format(np.mean(auc_list),np.std(auc_list)))\n",
    "print('mean precision={}±{}'.format(np.mean(precision_list),np.std(precision_list)))\n",
    "print('mean recall={}±{}'.format(np.mean(recall_list),np.std(recall_list)))\n",
    "print('mean f1_score={}±{}'.format(np.mean(f1_score_list),np.std(f1_score_list)))\n",
    "print('mean threshold={}±{}'.format(np.mean(thresh_list),np.std(thresh_list)))\n",
    "\n",
    "print('=========================================================')\n",
    "\n",
    "print('mean accuracy independent={}±{}'.format(np.mean(accuracy_list_t),np.std(accuracy_list_t)))\n",
    "print('mean auc independent={}±{}'.format(np.mean(auc_list_t),np.std(auc_list_t)))\n",
    "print('mean precision independent={}±{}'.format(np.mean(precision_list_t),np.std(precision_list_t)))\n",
    "print('mean recall independent={}±{}'.format(np.mean(recall_list_t),np.std(recall_list_t)))\n",
    "print('mean f1_score independent={}±{}'.format(np.mean(f1_score_list_t),np.std(f1_score_list_t)))\n",
    "print('mean threshold independent={}±{}'.format(np.mean(thresh_list_t),np.std(thresh_list_t)))\n",
    "'''\n",
    "print('GEO mean accuracy={}±{}'.format(np.mean(accuracy_list_geo),np.std(accuracy_list_geo)))\n",
    "print('GEO mean auc={}±{}'.format(np.mean(auc_list_geo),np.std(auc_list_geo)))\n",
    "print('GEO mean precision={}±{}'.format(np.mean(precision_list_geo),np.std(precision_list_geo)))\n",
    "print('GEO mean recall={}±{}'.format(np.mean(recall_list_geo),np.std(recall_list_geo)))\n",
    "print('GEO mean f1_score={}±{}'.format(np.mean(f1_score_list_geo),np.std(f1_score_list_geo)))\n",
    "print('GEO mean threshold={}±{}'.format(np.mean(thresh_list_geo),np.std(thresh_list_geo)))\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77f870e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "auc_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdc21ee8",
   "metadata": {},
   "outputs": [],
   "source": [
    "auc_list_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8d2c0ed",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2094a1c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a8e0609",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "cdcbef73",
   "metadata": {},
   "source": [
    "\n",
    "## Grad-Cam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7cb46b2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image, display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c13adbe9",
   "metadata": {},
   "outputs": [],
   "source": [
    "img_size = (49, 49)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0445f950",
   "metadata": {},
   "outputs": [],
   "source": [
    "last_conv_layer_name = \"conv2d_2\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7bec3079",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_img_array(img_path, size):\n",
    "    # `img` is a PIL image of size 49x49\n",
    "    img = keras.preprocessing.image.load_img(img_path, target_size=size)\n",
    "    # `array` is a float32 Numpy array of shape (49, 49, 1)\n",
    "    array = keras.preprocessing.image.img_to_array(img)\n",
    "    # We add a dimension to transform our array into a \"batch\"\n",
    "    # of size (1, 49, 49, 1)\n",
    "    array = np.expand_dims(array, axis=0)\n",
    "    return array\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9ee4f07",
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_gradcam_heatmap(img_array, model, last_conv_layer_name, pred_index=None):\n",
    "    # First, we create a model that maps the input image to the activations\n",
    "    # of the last conv layer as well as the output predictions\n",
    "    grad_model = tf.keras.models.Model(\n",
    "        [model.inputs], [model.get_layer(last_conv_layer_name).output, model.output]\n",
    "    )\n",
    "\n",
    "    # Then, we compute the gradient of the top predicted class for our input image\n",
    "    # with respect to the activations of the last conv layer\n",
    "    with tf.GradientTape() as tape:\n",
    "        last_conv_layer_output, preds = grad_model(img_array)\n",
    "        if pred_index is None:\n",
    "            pred_index = tf.argmax(preds[0])\n",
    "        class_channel = preds[:, pred_index]\n",
    "\n",
    "    # This is the gradient of the output neuron (top predicted or chosen)\n",
    "    # with regard to the output feature map of the last conv layer\n",
    "    grads = tape.gradient(class_channel, last_conv_layer_output)\n",
    "\n",
    "    # This is a vector where each entry is the mean intensity of the gradient\n",
    "    # over a specific feature map channel\n",
    "    pooled_grads = tf.reduce_mean(grads, axis=(0, 1, 2))\n",
    "\n",
    "    # We multiply each channel in the feature map array\n",
    "    # by \"how important this channel is\" with regard to the top predicted class\n",
    "    # then sum all the channels to obtain the heatmap class activation\n",
    "    last_conv_layer_output = last_conv_layer_output[0]\n",
    "    heatmap = last_conv_layer_output @ pooled_grads[..., tf.newaxis]\n",
    "    heatmap = tf.squeeze(heatmap)\n",
    "\n",
    "    # For visualization purpose, we will also normalize the heatmap between 0 & 1\n",
    "    heatmap = tf.maximum(heatmap, 0) / tf.math.reduce_max(heatmap)\n",
    "    return heatmap.numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5f90541",
   "metadata": {},
   "outputs": [],
   "source": [
    "img_array=X.reshape(X.shape[0],49,49,1)[0]\n",
    "img_array.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54c52f2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.cm as c\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.xticks([])  #Remove the values of the x-axis\n",
    "plt.yticks([])  #Remove the values of the y-axis\n",
    "plt.imshow(img_array,c.Spectral)\n",
    "plt.savefig('D:/Mydata/figs/pathway_block.svg', dpi=600, format='svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d609589",
   "metadata": {},
   "outputs": [],
   "source": [
    "img_array_t=X_test.reshape(X_test.shape[0],49,49,1)[0]\n",
    "img_array_t.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "311f1eaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 10))\n",
    "plt.imshow(img_array_t[8:14,0:7],c.PiYG)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d629b271",
   "metadata": {},
   "outputs": [],
   "source": [
    "img_array=img_array.reshape(1,49,49,1)\n",
    "model=best_model\n",
    "last_conv_layer_name='conv2d_2'\n",
    "pred_index=0\n",
    "model.layers[-1].activation=None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4964523c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate class activation heatmap\n",
    "heatmap = make_gradcam_heatmap(img_array, model, last_conv_layer_name,pred_index)\n",
    "\n",
    "# Display heatmap\n",
    "plt.imshow(heatmap,c.PiYG)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2313a592",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_model.predict(img_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "105381d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ad62e9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.columns=['label','Sam_Name']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d106a4d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_positive=df.loc[df.label==1.0]\n",
    "df_positive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ac38aaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_negative=df.loc[df.label==0.0]\n",
    "df_negative"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39290feb",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_positive = np.array([np.load(dir_name+s+\".npy\") for s in df_positive.index]) \n",
    "X_negative = np.array([np.load(dir_name+s+\".npy\") for s in df_negative.index])\n",
    "print(X_positive.shape)\n",
    "print(X_negative.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b6429fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "heatmap_negative=[]\n",
    "for i in range(X_negative.shape[0]):\n",
    "    img_array=X_negative[i]\n",
    "    img_array = np.expand_dims(img_array, axis=0)\n",
    "    model=best_model\n",
    "    model.layers[-1].activation=None\n",
    "    heatmap= make_gradcam_heatmap(img_array, model, last_conv_layer_name)\n",
    "    heatmap_negative.append(heatmap)\n",
    "    \n",
    "heatmap_negative=np.array(heatmap_negative)\n",
    "print(heatmap_negative.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37d5f67c",
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_mean_n=np.mean(heatmap_negative,axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c53424aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_mean_n.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f78bf274",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 10))\n",
    "plt.imshow(hm_mean_n,c.Reds)\n",
    "plt.xticks([])  \n",
    "plt.yticks([])  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f80efa1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "heatmap_positive=[]\n",
    "for i in range(X_positive.shape[0]):\n",
    "    img_array=X_positive[i]\n",
    "    img_array = np.expand_dims(img_array, axis=0)\n",
    "    model=best_model\n",
    "    model.layers[-1].activation=None\n",
    "    heatmap= make_gradcam_heatmap(img_array, model, last_conv_layer_name)\n",
    "    heatmap_positive.append(heatmap)\n",
    "    \n",
    "heatmap_positive=np.array(heatmap_positive)\n",
    "print(heatmap_positive.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61078472",
   "metadata": {},
   "outputs": [],
   "source": [
    "hm_mean_p=np.mean(heatmap_positive,axis=0)\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.imshow(hm_mean_p,c.Reds)\n",
    "plt.xticks([])  \n",
    "plt.yticks([])  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4eab9e0b",
   "metadata": {},
   "source": [
    "## Key Pathways Finding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36def9fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69f8ec37",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_hm_p=pd.DataFrame(heatmap_positive.reshape(heatmap_positive.shape[0],49))\n",
    "\n",
    "\n",
    "df_hm_p.columns=['Androgen_receptor_signaling_pathway','Cytokine-Jak-STAT_signaling_pathway','DCC-apoptotic_pathway','EGF-EGFR-PI3K_signaling_pathway','EGF-EGFR-PLCG-CAMK_signaling_pathway','EGF-EGFR-RAS-ERK_signaling_pathway','EGF-EGFR-RAS-RalGDS_signaling_pathway','EGF-EGFR-RAS-RASSF1_signaling_pathway','EGF-ERBB2-RAS-ERK_signaling_pathway','EGF-Jak-STAT_signaling_pathway','FGF-FGFR-PI3K_signaling_pathway','FGF-FGFR-RAS-ERK_signaling_pathway','FLT3LG-FLT3-RAS-ERK_signaling_pathway','FLT3LG-FLT3-RAS-PI3K_signaling_pathway','FLT3LG-FLT3-STAT5_signaling_pathway','Hedgehog_signaling_pathway','HGF-MET-PI3K_signaling_pathway','HGF-MET-RAS-ERK_signaling_pathway','HIF-1_signaling_pathway','IGF2-IGF1R-PI3K_signaling_pathway','IGF-IGF1R-RAS-ERK_signaling_pathway','Intrinsic_apoptotic_pathway','KEAP1-NRF2_signaling_pathway','KITLG-KIT-PI3K_signaling_pathway','KITLG-KIT-RAS-ERK_signaling_pathway','MDM2-p21-Cell_cycle_G1-S','NOTCH_signaling_pathway','p15-Cell_cycle_G1-S','p16-Cell_cycle_G1-S','p27-Cell_cycle_G1-S','PDGF-PDGFR-PI3K_signaling_pathway','PDGF-PDGFR-PLCG-CAMK_signaling_pathway','PDGF-PDGFR-RAS-ERK_signaling_pathway','PTEN-PIP3-AKT_signaling_pathway','Telomerase_activity','TGFA-EGFR-PI3K_signaling_pathway','TGFA-EGFR-PLCG-PKC_signaling_pathway','TGFA-EGFR-RAS-ERK_signaling_pathway','TGF-beta_signaling_pathway','WNT_signaling_pathway','z1','z2','z3','z4','z5','z6','z7','z8','z9']\n",
    "\n",
    "\n",
    "\n",
    "df_hm_p.drop('z1',axis=1,inplace=True)\n",
    "\n",
    "df_hm_p.drop('z2',axis=1,inplace=True)\n",
    "\n",
    "\n",
    "df_hm_p.drop('z3',axis=1,inplace=True)\n",
    "\n",
    "df_hm_p.drop('z4',axis=1,inplace=True)\n",
    "\n",
    "\n",
    "df_hm_p.drop('z5',axis=1,inplace=True)\n",
    "\n",
    "\n",
    "df_hm_p.drop('z6',axis=1,inplace=True)\n",
    "\n",
    "\n",
    "df_hm_p.drop('z7',axis=1,inplace=True)\n",
    "\n",
    "df_hm_p.drop('z8',axis=1,inplace=True)\n",
    "\n",
    "\n",
    "df_hm_p.drop('z9',axis=1,inplace=True)\n",
    "df_hm_p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14c69a99",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_hm_n=pd.DataFrame(heatmap_negative.reshape(heatmap_negative.shape[0],49))\n",
    "\n",
    "\n",
    "df_hm_n.columns=['Androgen_receptor_signaling_pathway','Cytokine-Jak-STAT_signaling_pathway','DCC-apoptotic_pathway','EGF-EGFR-PI3K_signaling_pathway','EGF-EGFR-PLCG-CAMK_signaling_pathway','EGF-EGFR-RAS-ERK_signaling_pathway','EGF-EGFR-RAS-RalGDS_signaling_pathway','EGF-EGFR-RAS-RASSF1_signaling_pathway','EGF-ERBB2-RAS-ERK_signaling_pathway','EGF-Jak-STAT_signaling_pathway','FGF-FGFR-PI3K_signaling_pathway','FGF-FGFR-RAS-ERK_signaling_pathway','FLT3LG-FLT3-RAS-ERK_signaling_pathway','FLT3LG-FLT3-RAS-PI3K_signaling_pathway','FLT3LG-FLT3-STAT5_signaling_pathway','Hedgehog_signaling_pathway','HGF-MET-PI3K_signaling_pathway','HGF-MET-RAS-ERK_signaling_pathway','HIF-1_signaling_pathway','IGF2-IGF1R-PI3K_signaling_pathway','IGF-IGF1R-RAS-ERK_signaling_pathway','Intrinsic_apoptotic_pathway','KEAP1-NRF2_signaling_pathway','KITLG-KIT-PI3K_signaling_pathway','KITLG-KIT-RAS-ERK_signaling_pathway','MDM2-p21-Cell_cycle_G1-S','NOTCH_signaling_pathway','p15-Cell_cycle_G1-S','p16-Cell_cycle_G1-S','p27-Cell_cycle_G1-S','PDGF-PDGFR-PI3K_signaling_pathway','PDGF-PDGFR-PLCG-CAMK_signaling_pathway','PDGF-PDGFR-RAS-ERK_signaling_pathway','PTEN-PIP3-AKT_signaling_pathway','Telomerase_activity','TGFA-EGFR-PI3K_signaling_pathway','TGFA-EGFR-PLCG-PKC_signaling_pathway','TGFA-EGFR-RAS-ERK_signaling_pathway','TGF-beta_signaling_pathway','WNT_signaling_pathway','z1','z2','z3','z4','z5','z6','z7','z8','z9']\n",
    "\n",
    "df_hm_n.drop('z1',axis=1,inplace=True)\n",
    "\n",
    "\n",
    "df_hm_n.drop('z2',axis=1,inplace=True)\n",
    "\n",
    "\n",
    "df_hm_n.drop('z3',axis=1,inplace=True)\n",
    "\n",
    "\n",
    "df_hm_n.drop('z4',axis=1,inplace=True)\n",
    "\n",
    "\n",
    "df_hm_n.drop('z5',axis=1,inplace=True)\n",
    "\n",
    "df_hm_n.drop('z6',axis=1,inplace=True)\n",
    "\n",
    "\n",
    "df_hm_n.drop('z7',axis=1,inplace=True)\n",
    "\n",
    "\n",
    "df_hm_n.drop('z8',axis=1,inplace=True)\n",
    "\n",
    "\n",
    "df_hm_n.drop('z9',axis=1,inplace=True)\n",
    "df_hm_n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d9a7828",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import ranksums\n",
    "\n",
    "pvalue_list=[]\n",
    "for col_p in df_hm_p.columns:\n",
    "    pvalue_list.append(ranksums(df_hm_p[col_p],df_hm_n[col_p]).pvalue)\n",
    "    print('Testing {}'.format(col_p))\n",
    "    print(ranksums(df_hm_p[col_p],df_hm_n[col_p]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a74df84",
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.stats.multitest import multipletests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4b350ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "multipletests(pvals=pvalue_list,alpha=0.05,method='bonferroni')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c664d425",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2a8f371",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "f86e2757",
   "metadata": {},
   "source": [
    "## Kaplan Meier Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2c83ef6",
   "metadata": {},
   "outputs": [],
   "source": [
    "pathways=pd.read_csv('D:/MyData/Pathway_Attention/data/All_Sub-pathways/allPathways.csv',sep='\\,',engine='python',header=None,index_col=0)\n",
    "pathways=pathways.fillna(0)\n",
    "pathways=pathways.drop(pathways.index[40])\n",
    "pathways"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13379676",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp=pd.read_hdf('data/SKCM/train_set.h5',key='gene_expression')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6670d8b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp.index.name=''\n",
    "exp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7dddeb02",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp.drop('geneName',axis=1,inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fa0e1f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp=exp.T\n",
    "exp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10c31d9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "tcga_clinical=pd.read_hdf('data/SKCM/train_set.h5',key='clinical')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1ec4123",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp.to_csv('tmp/exp_tcga_skcm.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1b40a24",
   "metadata": {},
   "outputs": [],
   "source": [
    "pathway_exp=pd.read_csv('tmp/res_tcga_skcm.csv',sep=',')\n",
    "pathway_exp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3339c7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "new_index=pathway_exp['Unnamed: 0']\n",
    "new_index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37828694",
   "metadata": {},
   "outputs": [],
   "source": [
    "pathway_exp=pathway_exp.iloc[:,1:578]\n",
    "pathway_exp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49d5e5f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "pathway_exp=pathway_exp.set_index(new_index)\n",
    "pathway_exp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47b2900e",
   "metadata": {},
   "outputs": [],
   "source": [
    "pathway_exp_new=pathway_exp.T.loc[[s for s in tcga_clinical[tcga_clinical['Pat_Died'].notna()].index]].T\n",
    "pathway_exp_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cb0cb32",
   "metadata": {},
   "outputs": [],
   "source": [
    "pathway_exp_median=pathway_exp_new.median(axis=1)\n",
    "pathway_exp_median.loc['Androgen_receptor_signaling_pathway']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4bde3ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "pathway_exp_mean=pathway_exp_new.mean(axis=1)\n",
    "pathway_exp_mean.loc['Androgen_receptor_signaling_pathway']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49766bab",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Testing Androgen_receptor_signaling_pathway\n",
    "RanksumsResult(statistic=-0.14786037739701, pvalue=0.8824529586738998)\n",
    "Testing Cytokine-Jak-STAT_signaling_pathway\n",
    "RanksumsResult(statistic=4.431752409589597, pvalue=9.3470310673971e-06)\n",
    "Testing DCC-apoptotic_pathway\n",
    "RanksumsResult(statistic=-1.665893585339646, pvalue=0.09573461149810288)\n",
    "Testing EGF-EGFR-PI3K_signaling_pathway\n",
    "RanksumsResult(statistic=-0.5746260156879879, pvalue=0.5655442534841102)\n",
    "Testing EGF-EGFR-PLCG-CAMK_signaling_pathway\n",
    "RanksumsResult(statistic=1.672851720746564, pvalue=0.09435649010598933)\n",
    "Testing EGF-EGFR-RAS-ERK_signaling_pathway\n",
    "RanksumsResult(statistic=0.9097762044545439, pvalue=0.36294054621992466)\n",
    "Testing EGF-EGFR-RAS-RalGDS_signaling_pathway\n",
    "RanksumsResult(statistic=-0.3739997781218488, pvalue=0.7084044812562552)\n",
    "Testing EGF-EGFR-RAS-RASSF1_signaling_pathway\n",
    "RanksumsResult(statistic=0.4215470367357893, pvalue=0.673355671682456)\n",
    "Testing EGF-ERBB2-RAS-ERK_signaling_pathway\n",
    "RanksumsResult(statistic=0.6546445728675462, pvalue=0.5126966191684176)\n",
    "Testing EGF-Jak-STAT_signaling_pathway\n",
    "RanksumsResult(statistic=-1.0686536295791742, pvalue=0.28522577334482835)\n",
    "Testing FGF-FGFR-PI3K_signaling_pathway\n",
    "RanksumsResult(statistic=1.6322625975395417, pvalue=0.10262417888697338)\n",
    "Testing FGF-FGFR-RAS-ERK_signaling_pathway\n",
    "RanksumsResult(statistic=2.7710774258051405, pvalue=0.0055871139815574745)\n",
    "Testing FLT3LG-FLT3-RAS-ERK_signaling_pathway\n",
    "RanksumsResult(statistic=-0.6070973142536058, pvalue=0.5437863324429668)\n",
    "Testing FLT3LG-FLT3-RAS-PI3K_signaling_pathway\n",
    "RanksumsResult(statistic=-2.4962310772318745, pvalue=0.012552080381874348)\n",
    "Testing FLT3LG-FLT3-STAT5_signaling_pathway\n",
    "RanksumsResult(statistic=-1.9871275032923656, pvalue=0.046908269114066564)\n",
    "Testing Hedgehog_signaling_pathway\n",
    "RanksumsResult(statistic=-1.2089760269520229, pvalue=0.22667205558953363)\n",
    "Testing HGF-MET-PI3K_signaling_pathway\n",
    "RanksumsResult(statistic=-3.6924505226045476, pvalue=0.00022210356044391689)\n",
    "Testing HGF-MET-RAS-ERK_signaling_pathway\n",
    "RanksumsResult(statistic=-1.3423402889179534, pvalue=0.17948568579748558)\n",
    "Testing HIF-1_signaling_pathway\n",
    "RanksumsResult(statistic=-0.7253856161712138, pvalue=0.46821549479831226)\n",
    "Testing IGF2-IGF1R-PI3K_signaling_pathway\n",
    "RanksumsResult(statistic=-0.7891685240679631, pvalue=0.43001351601890025)\n",
    "Testing IGF-IGF1R-RAS-ERK_signaling_pathway\n",
    "RanksumsResult(statistic=-4.177780467237087, pvalue=2.943675680933014e-05)\n",
    "Testing Intrinsic_apoptotic_pathway\n",
    "RanksumsResult(statistic=-1.1150411989586284, pvalue=0.26483276730645366)\n",
    "Testing KEAP1-NRF2_signaling_pathway\n",
    "RanksumsResult(statistic=-1.3516178027938444, pvalue=0.17649761218652582)\n",
    "Testing KITLG-KIT-PI3K_signaling_pathway\n",
    "RanksumsResult(statistic=-2.049750721954629, pvalue=0.04038876258568114)\n",
    "Testing KITLG-KIT-RAS-ERK_signaling_pathway\n",
    "RanksumsResult(statistic=0.13974255275560551, pvalue=0.8888634040158399)\n",
    "Testing MDM2-p21-Cell_cycle_G1-S\n",
    "RanksumsResult(statistic=-0.9352893676132437, pvalue=0.3496391932337082)\n",
    "Testing NOTCH_signaling_pathway\n",
    "RanksumsResult(statistic=-1.3087093011178492, pvalue=0.19063283988233726)\n",
    "Testing p15-Cell_cycle_G1-S\n",
    "RanksumsResult(statistic=-2.004522841809661, pvalue=0.04501408259253288)\n",
    "Testing p16-Cell_cycle_G1-S\n",
    "RanksumsResult(statistic=-1.2530442178625043, pvalue=0.21018961496983957)\n",
    "Testing p27-Cell_cycle_G1-S\n",
    "RanksumsResult(statistic=-2.6011829529528896, pvalue=0.009290289397190189)\n",
    "Testing PDGF-PDGFR-PI3K_signaling_pathway\n",
    "RanksumsResult(statistic=-0.8123623087576902, pvalue=0.41658376701889055)\n",
    "Testing PDGF-PDGFR-PLCG-CAMK_signaling_pathway\n",
    "RanksumsResult(statistic=1.162588457572569, pvalue=0.24499651607477113)\n",
    "Testing PDGF-PDGFR-RAS-ERK_signaling_pathway\n",
    "RanksumsResult(statistic=-2.31647924588649, pvalue=0.020532116329718415)\n",
    "Testing PTEN-PIP3-AKT_signaling_pathway\n",
    "RanksumsResult(statistic=-1.6357416652430008, pvalue=0.10189366543743018)\n",
    "Testing Telomerase_activity\n",
    "RanksumsResult(statistic=-1.7934594011331448, pvalue=0.07289949675453068)\n",
    "Testing TGFA-EGFR-PI3K_signaling_pathway\n",
    "RanksumsResult(statistic=-2.576249634411433, pvalue=0.009987850825986369)\n",
    "Testing TGFA-EGFR-PLCG-PKC_signaling_pathway\n",
    "RanksumsResult(statistic=-3.2384321873031405, pvalue=0.001201885892535581)\n",
    "Testing TGFA-EGFR-RAS-ERK_signaling_pathway\n",
    "RanksumsResult(statistic=-1.8815957829541077, pvalue=0.059890919240914685)\n",
    "Testing TGF-beta_signaling_pathway\n",
    "RanksumsResult(statistic=-0.9329699891442709, pvalue=0.3508354636604267)\n",
    "Testing WNT_signaling_pathway\n",
    "RanksumsResult(statistic=2.3895396676591303, pvalue=0.01686950203034994)\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f6467c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "pathway='IGF-IGF1R-RAS-ERK_signaling_pathway'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe696e8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "over_pat_list=[]\n",
    "lower_pat_list=[]\n",
    "for pat in pathway_exp.columns:\n",
    "    if pat!='geneName':\n",
    "        if pathway_exp[pat].loc[pathway]>=pathway_exp_median.loc[pathway]:\n",
    "            over_pat_list.append(pat)\n",
    "        if pathway_exp[pat].loc[pathway]<pathway_exp_median.loc[pathway]:\n",
    "            lower_pat_list.append(pat)\n",
    "print(len(over_pat_list))\n",
    "print(len(lower_pat_list))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9328f07e",
   "metadata": {},
   "outputs": [],
   "source": [
    "tcga_clinical_over=tcga_clinical[tcga_clinical.index.isin(over_pat_list)]\n",
    "tcga_clinical_lower=tcga_clinical[tcga_clinical.index.isin(lower_pat_list)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16930a1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "tcga_clinical_over=tcga_clinical_over[tcga_clinical_over['Pat_Overall_Survival_Months'].notna()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01797264",
   "metadata": {},
   "outputs": [],
   "source": [
    "tcga_clinical_lower=tcga_clinical_lower[tcga_clinical_lower['Pat_Overall_Survival_Months'].notna()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "350e9b2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(tcga_clinical_over)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0cdbae76",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(tcga_clinical_lower)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e755c7a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "from lifelines import KaplanMeierFitter\n",
    "timeline=[]\n",
    "for i in range(0,120):\n",
    "    timeline.append(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "582af2f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "from lifelines.statistics import logrank_test\n",
    "from lifelines.utils import restricted_mean_survival_time\n",
    "from lifelines.plotting import rmst_plot\n",
    "\n",
    "T_p=tcga_clinical_over['Pat_Overall_Survival_Months']\n",
    "E_p=tcga_clinical_over['Pat_Died']\n",
    "kmf_positive = KaplanMeierFitter()\n",
    "kmf_positive.fit(durations=T_p, event_observed=E_p,timeline=timeline,label='high expression')\n",
    "ax=kmf_positive.plot()\n",
    "\n",
    "T_n=tcga_clinical_lower['Pat_Overall_Survival_Months']\n",
    "E_n=tcga_clinical_lower['Pat_Died']\n",
    "kmf_negative=KaplanMeierFitter()\n",
    "kmf_negative.fit(durations=T_n, event_observed=E_n,timeline=timeline,label='low expression')\n",
    "ax=kmf_negative.plot(ax=ax)\n",
    "\n",
    "results = logrank_test(T_p, T_n, event_observed_A=E_p, event_observed_B=E_n)\n",
    "'''\n",
    "rmst_positive= restricted_mean_survival_time(kmf_positive, t=900)\n",
    "rmst_negative= restricted_mean_survival_time(kmf_negative, t=900)\n",
    "\n",
    "ax = plt.subplot(311)\n",
    "rmst_plot(kmf_negative, t=900, ax=ax)\n",
    "\n",
    "\n",
    "ax = plt.subplot(312)\n",
    "rmst_plot(kmf_positive, t=900, ax=ax)\n",
    "\n",
    "\n",
    "ax = plt.subplot(313)\n",
    "rmst_plot(kmf_negative, model2=kmf_positive, t=900, ax=ax)\n",
    "'''\n",
    "plt.legend(loc='lower left', bbox_to_anchor=(0.01,0.01),ncol=1,fancybox=False,shadow=False)\n",
    "plt.title('FLT3LG-FLT3-RAS-PI3K_signaling_pathway p-value: {}'.format(results.p_value))\n",
    "#plt.title('TGFA-EGFR-PLCG-PKC Signaling Pathway p-value = 0.1')\n",
    "ax.set_ylabel('Survival Rate')\n",
    "ax.set_xlabel('Timeline(months)')\n",
    "#ax.get_figure().savefig(\"FLT3LG-FLT3-RAS-PI3K_signaling_pathway.svg\",dpi=500,format='svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41a1281a",
   "metadata": {},
   "outputs": [],
   "source": [
    "km_list=[\n",
    "  0.6263281828023903,\n",
    "    1.2714239306329612e-05,\n",
    "    0.007211388960999662,\n",
    "    0.04485100190654674,\n",
    "    0.040436798213827736,\n",
    "    0.5109019953189557,\n",
    "    0.9615263409301946,\n",
    "    0.09030710787769194,\n",
    "    0.3917560260170282,\n",
    "    0.016958745259708023,\n",
    "    0.022028482481060985,\n",
    "    0.7422855498536335,\n",
    "    0.18962747499204072,\n",
    "    0.3680210492294468,\n",
    "    0.0005649350473243639,\n",
    "    0.13882123394186013,\n",
    "    0.045679067370708974,\n",
    "    0.4424199375183385,\n",
    "    0.30712483615717756,\n",
    "    0.6453481450542624,\n",
    "    0.10432074129391562,\n",
    "    0.505477325178186,\n",
    "    0.0007550269096438157,\n",
    "    0.00022797808712110318,\n",
    "    0.11654609027891244,\n",
    "    0.025863734633726664,\n",
    "    0.0004707197071112644,\n",
    "    0.0014393233463285123,\n",
    "    0.08829142376075509,\n",
    "    0.026610000015394552,\n",
    "    0.1383310708761253,\n",
    "    0.24390443553712657,\n",
    "    0.5133981996044519,\n",
    "    0.18021362756990777,\n",
    "    0.3316838212080421,\n",
    "    0.3197032087390901,\n",
    "    0.10973166123963331,\n",
    "    0.21138495173337732,\n",
    "    0.6593409669277395,\n",
    "    0.05863179875693078\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb0ebae2",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(km_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d3445b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "bonferroni_list=[1.00000000e+00, 3.73881243e-04, 1.00000000e+00, 1.00000000e+00,\n",
    "        1.00000000e+00, 1.00000000e+00, 1.00000000e+00, 1.00000000e+00,\n",
    "        1.00000000e+00, 1.00000000e+00, 1.00000000e+00, 2.23484559e-01,\n",
    "        1.00000000e+00, 5.02083215e-01, 1.00000000e+00, 1.00000000e+00,\n",
    "        8.88414242e-03, 1.00000000e+00, 1.00000000e+00, 1.00000000e+00,\n",
    "        1.17747027e-03, 1.00000000e+00, 1.00000000e+00, 1.00000000e+00,\n",
    "        1.00000000e+00, 1.00000000e+00, 1.00000000e+00, 1.00000000e+00,\n",
    "        1.00000000e+00, 3.71611576e-01, 1.00000000e+00, 1.00000000e+00,\n",
    "        8.21284653e-01, 1.00000000e+00, 1.00000000e+00, 3.99514033e-01,\n",
    "        5.2754357e-02, 1.00000000e+00, 1.00000000e+00, 6.74780081e-01,1,1,1,1,1,1,1,1,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26baa23a",
   "metadata": {},
   "outputs": [],
   "source": [
    "salient_array=np.array(bonferroni_list).reshape(7,7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77225d18",
   "metadata": {},
   "outputs": [],
   "source": [
    "salient_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de0970ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 14))\n",
    "plt.xticks([])  \n",
    "plt.yticks([])  \n",
    "plt.imshow(salient_array,c.afmhot)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68ac12f3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}