DeepLearning_CyTOF.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center> <h1> Diagnose Latent Cytomegalovirus Using CyTOF Data and Deep Learning <h1> </center>\n",
    "\n",
    "<center>Zicheng Hu, Ph.D.</center>\n",
    "<center>Research Scientist</center>\n",
    "<center>ImmPort Team</center>\n",
    "<center>The Unversity of California, San Francisco</center>\n",
    "\n",
    "![alt text](Data/header.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction\n",
    "A deep neural network (a.k.a. deep learning) is an artificial neural network with multiple layers between the input and output layers. It was proven to be highly effective for a variety of predictive tasks. In health care, deep learning is quickly gaining popularity and has been implemented for applications such as image-based diagnosis and personalized drug recommendations. In this tutorial, we will build a tailored deep-learning model for CyTOF data to diagnosis latent Cytomegalovirus infection using Keras and TensorFlow. To run this tutorial, download the [github repository](https://github.com/hzc363/DeepLearningCyTOF) and run the [jupyter notebook](https://github.com/hzc363/DeepLearningCyTOF/blob/master/DeepLearning_CyTOF.ipynb). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: Import Functions\n",
    "Before we start, we first import functions that we will use in this tutorial from different libraries. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### Step 1: import functions #####\n",
    "from keras.layers import Dense, Flatten, BatchNormalization, Activation, Conv2D, AveragePooling2D, Input\n",
    "from keras.models import load_model, Model\n",
    "from keras.optimizers import Adam\n",
    "from keras.callbacks import ModelCheckpoint, EarlyStopping\n",
    "from keras import backend as K\n",
    "import pickle\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from numpy.random import seed; seed(111)\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from tensorflow import set_random_seed; set_random_seed(111)\n",
    "from sklearn.metrics import roc_curve, auc\n",
    "from sklearn.externals.six import StringIO  \n",
    "from sklearn.tree import export_graphviz, DecisionTreeRegressor\n",
    "from scipy.stats import ttest_ind\n",
    "from IPython.display import Image  \n",
    "import pydotplus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: Load data\n",
    "\n",
    "We load the data, which are stored in the \"allData.obj\" file. The data includes three parts, meta-data, CyTOF data, and marker names. \n",
    "\n",
    "* The **CyTOF data** contains the single-cell profile of 27 markers. For the convenience of this tutorial, we already downloaded the fcs files from ImmPort and preprocessed the data into Numpy arrays. See [an example](https://github.com/hzc363/DeepLearningCyTOF/tree/master/FCS_to_Array) for the preprocessing of the FCS files. The dimension of the Numpy array is 472 samples x 10000 cells x 27 markers.\n",
    "* The **metadata** contains the sample level information, including the study accession number for each sample and the ground truth of CMV infection. It is stored as a pandas data frame.\n",
    "* The **marker names** contain the name of the 27 markers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### Step 2: load data #####\n",
    "\n",
    "#Download data\n",
    "tutorial_files = ! ls Data\n",
    "if \"allData.obj\" not in tutorial_files:\n",
    "    print(\"Downloading Data:\")\n",
    "    ! wget https://figshare.com/ndownloader/files/38918480 -P ./Data\n",
    "        \n",
    "#load data\n",
    "allData = pickle.load( open( \"Data/allData.obj\", \"rb\" ) )\n",
    "metaData = allData[\"metaData\"]\n",
    "cytoData = allData[\"cytoData\"]\n",
    "markerNames = allData[\"markerNames\"]\n",
    "\n",
    "# inspect the data\n",
    "print(\"\\nFirst 5 rows of metaData: \")\n",
    "print(metaData.head(),\"\\n\")\n",
    "\n",
    "print(\"Dimensions of cytoData: \",cytoData.shape,\"\\n\")\n",
    "print(\"Names of the 27 makers: \\n\",markerNames.values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3: Split data into training, validation and testing sets\n",
    "Now, lets split the data into training, validation, and testing sets. The training data is used to train the deep learning model. The validation dataset is used to select the best parameters for the model and to avoid overfitting. The test dataset is used to evaluate the performance of the final model.\n",
    "\n",
    "The CyTOF dataset contains samples from 9 studies available on ImmPort. We will use samples from the study SDY515 as a validation set, samples from the study SDY519 as a testing set, and the rest of the samples as a training set. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### Step 3: split train, validation and test######\n",
    "y = metaData.CMV_Ab.values\n",
    "x = cytoData\n",
    "\n",
    "train_id = (metaData.study_accession.isin([\"SDY515\",\"SDY519\"])==False)\n",
    "valid_id = metaData.study_accession==\"SDY515\"\n",
    "test_id = metaData.study_accession ==\"SDY519\"\n",
    "\n",
    "x_train = x[train_id]; y_train = y[train_id]\n",
    "x_valid = x[valid_id]; y_valid = y[valid_id]\n",
    "x_test = x[test_id]; y_test = y[test_id]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 4: Define the deep learning model\n",
    "We will use a customized convolution neural network (CNN) to analyze the CyTOF data. For each sample, the CyTOF data is a matrix with rows as cells and columns as markers. It is crucial to notice that the CyTOF data is an unordered collection of cells (rows). For example, both matrix 1 and matrix 2 profiles the same sample in Figure 1A, even though they have different orders of rows. \n",
    "\n",
    "![alt text](Data/Figure1.png)\n",
    "\n",
    "\n",
    "Based on the characteristics of the CyTOF data, we design a CNN model that is invariant to the permutation of rows. The model contains six layers: input layer, first and second convolution layer, pooling layer, dense layer, and output layer. \n",
    "\n",
    "* The **input layer** receives the CyTOF data matrix. \n",
    "\n",
    "* The **first convolution layer** uses three filters to scan each row of the CyTOF data. This layer extracts relevant information from the cell marker profile of each cell. \n",
    "\n",
    "* The **second convolution layer** uses three filters to scan each row of the first layer's output. Each filter combines information from the first layer for each cell. \n",
    "\n",
    "* The **pooling layers** averages the outputs of the second convolution layer. The purpose is to aggregate the cell level information into sample-level information. \n",
    "\n",
    "* The **dense layer** further extracts information from the pooling layer. \n",
    "\n",
    "* The **output layer** uses logistic regression to report the probability of CMV infection for each sample. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### Step 4: define model #####\n",
    "\n",
    "# input\n",
    "model_input = Input(shape=x_train[0].shape)\n",
    "\n",
    "# first convolution layer\n",
    "model_output = Conv2D(3, kernel_size=(1, x_train.shape[2]),\n",
    "                 activation=None)(model_input)\n",
    "model_output = BatchNormalization()(model_output)\n",
    "model_output = Activation(\"relu\")(model_output)\n",
    "\n",
    "# sceond convolution layer\n",
    "model_output = Conv2D(3, (1, 1), activation=None)(model_output)\n",
    "model_output = BatchNormalization()(model_output)\n",
    "model_output = Activation(\"relu\")(model_output)\n",
    "\n",
    "# pooling layer\n",
    "model_output = AveragePooling2D(pool_size=(x_train.shape[1], 1))(model_output)\n",
    "model_output = Flatten()(model_output)\n",
    "\n",
    "# Dense layer\n",
    "model_output = Dense(3, activation=None)(model_output)\n",
    "model_output = BatchNormalization()(model_output)\n",
    "model_output = Activation(\"relu\")(model_output)\n",
    "\n",
    "# output layer\n",
    "model_output = Dense(1, activation=None)(model_output)\n",
    "model_output = BatchNormalization()(model_output)\n",
    "model_output = Activation(\"sigmoid\")(model_output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 5: Fit the model\n",
    "In this step, we will use the training data to fit the model. We will use the Adam algorithm, which is an extension of the gradient descent method to train our model. Adam algorithm will search the model space step by step (epochs) until the optimal model is identified. At each step, we will use validation data to evaluate the performance of the model. The best model will be saved. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### Step 5: Fit model #####\n",
    "\n",
    "# specify input and output\n",
    "model = Model(inputs=[model_input],\n",
    "              outputs=model_output)\n",
    "\n",
    "# define loss function and optimizer\n",
    "model.compile(loss='binary_crossentropy',\n",
    "              optimizer=Adam(lr=0.0001),\n",
    "              metrics=['accuracy'])\n",
    "\n",
    "# save the best performing model\n",
    "checkpointer = ModelCheckpoint(filepath='Result/saved_weights.hdf5', \n",
    "                               monitor='val_loss', verbose=0, \n",
    "                               save_best_only=True)\n",
    "\n",
    "# model training\n",
    "model.fit([x_train], y_train,\n",
    "          batch_size=60,\n",
    "          epochs=500, \n",
    "          verbose=1,\n",
    "          callbacks=[checkpointer],\n",
    "          validation_data=([x_valid], y_valid))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 6: Plot the training history\n",
    "We can view the training history of the model by plotting the performance (value of the loss function) for training and validation data in each epoch. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "##### Step 6: plot train and validation loss #####\n",
    "plt.plot(model.history.history['loss'])\n",
    "plt.plot(model.history.history['val_loss'])\n",
    "plt.title('model train vs validation loss')\n",
    "plt.ylabel('loss')\n",
    "plt.xlabel('epoch')\n",
    "plt.legend(['train', 'validation'], loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 7: Evaluate the performance using test data\n",
    "We load the final model from a save file (Final_weights.hdf5) for the following analysis steps. We will use the test data, which has not been touched so far, to evaluate the performance of the final model. We will draw a Receiver Operator Characteristic(ROC) Curve and use Area Under the Curve (AUC) to measure performance. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### Step 7: test the final model #####\n",
    "\n",
    "# load final model\n",
    "final_model = load_model('Data/Final_weights.hdf5')\n",
    "\n",
    "# generate ROC and AUC\n",
    "y_scores = final_model.predict([x_test])\n",
    "fpr, tpr, _ = roc_curve(y_test, y_scores)\n",
    "roc_auc = auc(fpr, tpr)\n",
    "\n",
    "# plot ROC curve\n",
    "plt.plot(fpr, tpr)\n",
    "plt.plot([0, 1], [0, 1], 'k--')\n",
    "plt.xlabel('False Positive Rate')\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.title('AUC = {0:.2f}'.format(roc_auc))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Step 8: Interpret the deep learning model.\n",
    "We will use a permutation-based method to interpret the deep CNN model (Fig. 2). For each cell in cytometry data, we up-sampled the cell by copying it to replace other randomly chosen cells within the sample. We then applied the deep CNN model on both the original data and the permuted data. The difference in the model output (ΔΥ) quantifies the impact of the cell on the output of the deep learning model. Finally, we build a decision tree model to idenitfy cell subsets that have the highest ΔΥ. \n",
    "\n",
    "![alt text](Data/Figure2.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### Step 8: Interpret the deep learning model. #####\n",
    "\n",
    "# warning: may take a long time (around 30 mins) to run\n",
    "\n",
    "# Calculate the impact of each cell on the model output\n",
    "dY = np.zeros([x_test.shape[0],x_test.shape[1]])\n",
    "s1 = np.random.randint(0,(x_test.shape[1]-1),int(x_test.shape[1]*0.05))\n",
    "\n",
    "for i in range(x_test.shape[0]):\n",
    "    pred_i = final_model.predict([x_test[[i],:,:,:]])\n",
    "    for j in range(x_test.shape[1]):\n",
    "        t1 = x_test[[i],:,:,:].copy()\n",
    "        t1[:,s1,:,:] = t1[:,j,:,:]\n",
    "        pred_j = final_model.predict([t1])\n",
    "        dY[i,j] = pred_j-pred_i\n",
    "\n",
    "# reformat dY\n",
    "x_test2 = x_test.reshape((x_test.shape[0]*x_test.shape[1],27))\n",
    "dY = dY.reshape([x_test.shape[0]*x_test.shape[1]])\n",
    "\n",
    "# Build decision tree to identify cell subset with high dY\n",
    "regr_1 = DecisionTreeRegressor(max_depth=4)\n",
    "regr_1.fit(x_test2, dY)\n",
    "\n",
    "# Plot the decision tree\n",
    "dot_data = StringIO()\n",
    "export_graphviz(regr_1, out_file=dot_data, \n",
    "                feature_names= markerNames,\n",
    "                filled=True, rounded=True,\n",
    "                special_characters=True)\n",
    "graph = pydotplus.graph_from_dot_data(dot_data.getvalue())  \n",
    "Image(graph.create_png())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conclusion\n",
    "In this tutorial, we build a deep convolutional neural network (CNN) to analyze CyTOF data. The deep CNN model is able to diagnose latent CMV infection with high accuracy. In addition, we were able to interpret the deep learning model using a permutation-based method. We discovered that a CD3+ CD8+ CD27- CD94+ population that have the highest impact on the deep CNN model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": 1
}