Workshop.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 0 - The Kaplan-Meyer Estimator\n",
    "\n",
    "Lets first familiarize ourselves with the Kaplan–Meier estimator. You might remember that its formula reads\n",
    "\n",
    "\n",
    "\\begin{equation*}\n",
    "\\widehat S(t) = \\prod_{i:\\ t_i\\le t} \\left( 1 - \\frac{d_i}{n_i} \\right)\n",
    "\\end{equation*}\n",
    "\n",
    "\n",
    "Where $\\widehat S(t)$ is the estimator for time $t$, $d_i$ the number of deaths at time $t_i$ and $n_i$ the number of individuals that are **known to have survived** until time $t_i$.\n",
    "\n",
    "With this in mind, use the following data to calculate $\\widehat S(t)$ and sketch a plot for the estimator.\n",
    "\n",
    "|Subject|Survival Time|Status|\n",
    "|---|---|---|\n",
    "|Fly 1|1 day|Dead|\n",
    "|Fly 2|2 days|Alive|\n",
    "|Fly 3|3 days|Alive|\n",
    "|Fly 4|4 days|Dead|\n",
    "|Fly 5|5 days|Dead|\n",
    "|Fly 6|6 days|Alive|\n",
    "|Fly 7|7 days|Dead|\n",
    "|Fly 8|8 days|Dead|\n",
    "|Fly 9|9 days|Alive|\n",
    "|Fly 10|10 days|Alive|\n",
    "\n",
    "\n",
    "\n",
    "---\n",
    "\n",
    "# Part 1 - Importing necessary packages\n",
    "\n",
    "This will import the packages we will use during the workshop. If the installation was done correctly, there should be no errors here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "# Numerical library\n",
    "import numpy as np\n",
    "\n",
    "# Data manipulation\n",
    "import pandas as pd\n",
    "from patsy import dmatrix\n",
    "\n",
    "# Ploting\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Survival analysis\n",
    "import lifelines\n",
    "\n",
    "from IPython.display import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 2 - A look at the clinical data\n",
    "\n",
    "The data is originally available at http://www.cbioportal.org/study?id=brca_tcga_pub, but it can also be found together with this notebook at https://github.com/gjeuken/survival-workshop.\n",
    "\n",
    "We will use here clinical data on the patients as well as gene expression data. This is data from real experiments (!) so you might even discover something that went unnoticed until now.\n",
    "\n",
    "Firstly, let us load the clinical data and see what it looks like. Make a note of all the variables (columns) available in the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "clinical = pd.read_csv('data/brca_tcga_pub_clinical_data.tsv', sep='\\t')\n",
    "\n",
    "with pd.option_context('display.max_columns', None):\n",
    "    display(clinical)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make the later survival analysis easier, we create two new columns, one for the survival in days, and one for the event (death) being present."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clinical['Time'] = clinical['Overall Survival (Months)']*30\n",
    "clinical['Event'] = (clinical['Overall Survival Status'] == 'DECEASED')*1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let us plot how the survival data looks like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "matplotlib.rcParams['figure.figsize'] = [15, 50]\n",
    "\n",
    "data_sorted = clinical[['Time','Event']].sort_values(by = 'Time').dropna().reset_index(drop=True)\n",
    "print(data_sorted.shape)\n",
    "status_slice = data_sorted['Event'] == 1\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.barh(data_sorted.loc[~status_slice].index, data_sorted.loc[~status_slice,'Time'], height = 1, color = 'b')\n",
    "ax.barh(data_sorted.loc[status_slice].index, data_sorted.loc[status_slice,'Time'], height = 1, color = 'r')\n",
    "\n",
    "ax.set_xlabel('Days alive')\n",
    "ax.set_ylabel('Patients')\n",
    "ax.legend(['Alive', 'Dead'])\n",
    "ax.xaxis.tick_top()\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Exercises\n",
    "2.1 Describe the plot and what inferences you are able to make from it\n",
    "\n",
    "---\n",
    "# Part 3 - Survival analysis using the clinical data\n",
    "\n",
    "\n",
    "\n",
    "Ok, now that we've seen the data, let's play around with it.\n",
    "\n",
    "How does the survival curve look like in general? We can use the survival package __lifelines__ to figure this out, and generate a *Kaplan Meier plot*. \n",
    "\n",
    "We will also compare the Kaplan Meier estimate for $S(t)$ with a naive estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "kmf = lifelines.KaplanMeierFitter()\n",
    "kmf_naive = lifelines.KaplanMeierFitter()\n",
    "\n",
    "clinical = clinical.dropna()\n",
    "\n",
    "kmf.fit(clinical['Time'], clinical['Event'])\n",
    "kmf_naive.fit(clinical['Time'], np.repeat(1, len(clinical['Time'])))\n",
    "\n",
    "matplotlib.rcParams['figure.figsize'] = [15, 10]\n",
    "ax = kmf.plot(ci_show=False)\n",
    "kmf_naive.plot(ax=ax, ci_show=False)\n",
    "plt.legend(['Kaplan-Meier estimate','Naive estimate'])\n",
    "plt.xlabel('Days')\n",
    "plt.ylabel('$S(t)$')\n",
    "plt.ylim([0,1])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Exercises\n",
    "\n",
    "3.1 Why are both estimates so different? Can you relate your answer to the first plot and the data?\n",
    "\n",
    "3.2 Compare the Kaplan Meier plot with the first one, what additional insights are available on this latter plot?\n",
    "\n",
    "3.3 What are the common features of the two plots?\n",
    "\n",
    "---\n",
    "\n",
    "Now we can start to play around with clinical variables that might influence in the survival curve.\n",
    "\n",
    "You should play around with the groupings and see if you can find some useful insight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Define groups here\n",
    "group_1 = clinical.loc[clinical['Metastasis-Coded'] == 'Positive']\n",
    "group_2 = clinical.loc[clinical['Metastasis-Coded'] == 'Negative']\n",
    "###\n",
    "\n",
    "kmf.fit(group_1['Time'], group_1['Event'], label='Positive')\n",
    "ax = kmf.plot(ci_show=False)\n",
    "\n",
    "kmf.fit(group_2['Time'], group_2['Event'], label='Negative')\n",
    "kmf.plot(ax=ax, ci_show=False)\n",
    "\n",
    "plt.title('Metastasis')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Exercises\n",
    "3.4 Make at least 3 plots using relevant separation criteria, and label them accordingly. Why did you choose these criteria? Are the results as expected?\n",
    "\n",
    "3.5 Make a plot that includes the survival curve for 3 different age groups. Do the results make sense?\n",
    "\n",
    "3.6 (advanced) The [**logrank_test** function](http://lifelines.readthedocs.io/en/latest/Examples.html?compare-using-a-hypothesis-test#compare-using-a-hypothesis-test) of the **lifelines** package performs a statistical test on the two groups to see if their event (death) generation process is the same. Use this function to obtain a significance statistic for the separations above.\n",
    "\n",
    "---\n",
    "\n",
    "# Part 4 - Cox Proportional Hazards model\n",
    "\n",
    "We will now introduce the Cox Proportional Hazards (CPH) model. To do so, lets first define the hazard function:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\lambda(t) = \\lim_{dt \\rightarrow 0} \\frac{\\Pr(t \\leq T < t+dt)}{dt\\cdot S(t)} = - \\frac{S'(t)}{S(t)}\n",
    "\\end{equation*}\n",
    "\n",
    "To understand the model, we first need to understand its assumptions. The main assumption here (proportional hazards condition) states that covariates are multiplicatively related to the hazard, or in other words, a change in the covariate will have an effect on the hazard rate that is proportional to the magnitude of that change.\n",
    "\n",
    "Now we can in introduce the model. The formula for CPH reads:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\lambda(t|X_i) = \\lambda_0(t)\\exp(\\beta_1X_{i1} + \\cdots + \\beta_pX_{ip}) = \\lambda_0(t)\\exp(X_i \\cdot \\beta)\n",
    "\\end{equation*}\n",
    "\n",
    "Where $X_i = \\{X_{i,1}, X_{i,2}, \\dots, X_{i,p}\\}$ are the realized values of the covariates for subject $i$, $\\lambda_0(t)$ is the *baseline* hazard rate, and $\\beta_i$ are the coefficients of proportionality.\n",
    "With a careful read, you can see how this equation relates to the assumption discussed. If you want to read more about the model, and especially how the coefficients are calculated, a great place to start is the [Wikipedia page](https://en.wikipedia.org/wiki/Proportional_hazards_model).\n",
    "\n",
    "The next step is to use the Cox Proportional Hazards regression in the data. \n",
    "First, let's check the main assumption of the model, the proportionality of hazards, by plotting the hazard functions on a log-log scale and seeing whether they \"look\" parallel or not. We will do this for a few covariates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Define groups here\n",
    "group_1 = clinical.loc[clinical['Metastasis-Coded'] == 'Positive']\n",
    "group_2 = clinical.loc[clinical['Metastasis-Coded'] == 'Negative']\n",
    "###\n",
    "\n",
    "kmf.fit(group_1['Time'], group_1['Event'], label='Positive')\n",
    "ax = kmf.plot_loglogs()\n",
    "\n",
    "kmf.fit(group_2['Time'], group_2['Event'], label='Negative')\n",
    "kmf.plot_loglogs(ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Exercises\n",
    "\n",
    "4.1 Can you explain what the baseline hazard means?\n",
    "\n",
    "4.2 Why can we test our hypothesis by looking at the plot above?\n",
    "\n",
    "4.3 In your opinion, does the proportional hazards assumption hold for the Metastasis case? What about other variables you used in the previous exercise? Use the plots to support your claims.\n",
    "\n",
    "---\n",
    "\n",
    "Now that we checked the proportional hazards assumption, we can perform the CPH regression on the categorical data where this assumption holds.\n",
    "\n",
    "Note that here we will use the [patsy package](https://patsy.readthedocs.io/en/latest/) to format our categorical data in a way that is suited for the regression. If you are not familiar with this type of representation, please take a moment to familiarize yourself in the [relevant help page](https://patsy.readthedocs.io/en/latest/categorical-coding.html).\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cph = lifelines.CoxPHFitter()\n",
    "\n",
    "formula = 'Time + Event + C(Q(\"Metastasis\"), levels = l) + Q(\"Diagnosis Age\")'\n",
    "\n",
    "l = ['M0', 'M1']\n",
    "\n",
    "\n",
    "cox_input = dmatrix(formula, data = clinical, return_type = 'dataframe')\n",
    "cox_input = cox_input.drop(['Intercept'], axis = 1)\n",
    "\n",
    "cph.fit(cox_input, duration_col='Time', event_col='Event')\n",
    "\n",
    "cph.plot()\n",
    "cph.summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Exercises\n",
    "4.4 For the regression performed above (Metastasis + Age), compare both coefficients. What can you learn from them?\n",
    "\n",
    "4.5 Perform the CPH regression using a few variables where the proportional hazard assumption hold.\n",
    "\n",
    "4.6 Try adding more than one variable in the same regression. What happens then? Explain the results.\n",
    "\n",
    "4.7 (advanced) If there are variables where the proportional hazard assumption does not hold, but hold predictive power, you might still use it for <a href=\"https://en.wikipedia.org/wiki/Blocking_(statistics)\">blocking</a>. How do the results change then? [Help](http://lifelines.readthedocs.io/en/latest/Survival%20Regression.html#stratification)\n",
    "\n",
    "4.8 Can you find a way to plot the survival curves predicted by the regression for each level of a categorical variable? (Remember the CPH assumptions)\n",
    "\n",
    "4.9 What are the difference in the results compared to the Kaplan-Mayer statistics when using the same covariates?\n",
    "\n",
    "---\n",
    "\n",
    "# Part 5 - Incorporating gene expression data\n",
    "\n",
    "It's time to look at the gene expression data! This data is generated using RNASeq on tumor tissue. \n",
    "\n",
    "First, we load the expression data for the same samples. This may take some time, be patient.\n",
    "\n",
    "We then take a look at how the data looks like. Note that the gene expression has been median normalized, so we get negative gene expression here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "expression_raw  = pd.read_csv('data/data_expression_median.txt', sep='\\t')\n",
    "expression = expression_raw.set_index('Hugo_Symbol').iloc[:,1:].T\n",
    "expression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we merge the clinical and expression data in one big data structure, and display the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "data = clinical.merge(expression, how='inner', left_on='Sample ID', right_index=True)\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Exercises\n",
    "\n",
    "5.1 How many genes are there in the data? Is there any gene you think might be relevant for the survival time in breast cancer?\n",
    "\n",
    "5.2 Can you think of natural ways to incorporate this data on both the Kaplan–Meier estimator and the Cox Proportional Hazards models?\n",
    "\n",
    "----\n",
    "\n",
    "Let's incorporate the data in our study. We will start with the Kaplan–Meier estimator, and to do so, we will have to separate our subjects in groups according to gene expression levels.\n",
    "We will next plot the Kaplan–Meier estimates, together with the results of the statistical test.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Define groups here\n",
    "treshold = 0\n",
    "\n",
    "group_1 = data.loc[data['BRCA1'] > treshold]\n",
    "group_2 = data.loc[data['BRCA1'] <= treshold]\n",
    "###\n",
    "\n",
    "kmf.fit(group_1['Time'], group_1['Event'], label='Positive')\n",
    "ax = kmf.plot(ci_show=False)\n",
    "\n",
    "kmf.fit(group_2['Time'], group_2['Event'], label='Negative')\n",
    "kmf.plot(ax=ax, ci_show=False)\n",
    "\n",
    "plt.title('BRCA1')\n",
    "\n",
    "from lifelines.statistics import logrank_test\n",
    "\n",
    "test = logrank_test(group_1['Time'], group_2['Time'], event_observed_A=group_1['Event'], event_observed_B=group_2['Event'])\n",
    "test.summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Exercises\n",
    "\n",
    "5.3 What criteria would you use to separate the groups using the gene expression information? Can you come up with a few?\n",
    "\n",
    "5.4 Perform the calculation with some genes of your choice, and using a few separation criteria. Do you get interesting results? Are they statistically significant?\n",
    "\n",
    "---\n",
    "\n",
    "We can also do the same for the CPH model. Let's incorporate the gene expression data into our Cox regression.\n",
    "\n",
    "Notice that here there is no need to separate the subjects into groups."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "formula = 'Time + Event + BRCA1'\n",
    "\n",
    "cox_input = dmatrix(formula, data = data, return_type = 'dataframe')\n",
    "cox_input = cox_input.drop(['Intercept'], axis = 1)\n",
    "\n",
    "cph.fit(cox_input, duration_col='Time', event_col='Event')\n",
    "\n",
    "cph.plot()\n",
    "cph.summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Exercises\n",
    "5.5 Perform the CPH regression with some genes of your choice. Did you find any gene that has a significant influence in the survival rate for breast cancer?\n",
    "\n",
    "5.6 How do the results differ from the test performed previously?\n",
    "\n",
    "5.7 (advanced) You may want to write a code that automatically tests all the genes and outputs all the relevant genes and their regression coefficients in the end.\n",
    "\n",
    "5.8 Now try multiple combinations of the genes in the same regression, start with the genes found in the previous exercise. Again, do the results change? Why?\n",
    "\n",
    "\n",
    "\n",
    "# The final exercise\n",
    "\n",
    "What is your opinion of survival analysis in general? How do you compare both models studied here? What are their merits and disadvantages? Where would you use one instead of the other?\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "----\n",
    "**Thank you for participating!** Any feedback on the content and delivery of the workshop is welcomed at gustavo.jeuken@scilifelab.se.\n"
   ]
  }
 ],
 "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.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}