Workshop_Notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "**Summary:** An introductory tutorial to get things started with single cell RNA-seq analysis. This tutorial will describe the basic steps in a single cell RNA-seq analysis workflow including cell and gene quality control, normalization, batch correction, selection of highly variable genes, and visualization. Once described, we will leverage these techniques to identify and characterize T cell subsets and determine differences in these populations that result from virus infection.\n",
    "\n",
    "Data required for workshop can be found at [Andersen Lab's Google Cloud](https://storage.cloud.google.com/andersen-lab_temp/dengue_all.hdf5). Download and move into the /data folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:01:21.873108Z",
     "start_time": "2020-02-24T19:01:15.653575Z"
    },
    "pycharm": {
     "is_executing": false
    }
   },
   "outputs": [],
   "source": [
    "import scanpy as sc\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import wget\n",
    "\n",
    "from rpy2.robjects import pandas2ri\n",
    "\n",
    "pandas2ri.activate()\n",
    "\n",
    "%load_ext rpy2.ipython\n",
    "\n",
    "sc.settings.set_figure_params( dpi=80 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Dataset Selection\n",
    "One benefit of the increased popularity of single cell RNA-seq is the release of a plentiful numbers of organism and cell type specific datasets. For this workshop, we'll be using a dataset provided by Zanini et al. from their paper [Virus-inclusive single-cell RNA sequencing reveals the molecular signature of progression to severe dengue](https://doi.org/10.1073/pnas.1813819115).\n",
    "\n",
    "In brief, blood samples were collected from four healthy controls and six individuals who presenting to a Colombian hospital and met criteria for dengue virus infection. Two of the dengue-infected patients experienced an uncomplicated disease course characterized by flu-like symptoms while the other four experienced severe dengue manifestions.  Peripheral blood mononuclear cells were extracted and FACS was performed to enrich various broad immune cell types and sort single cells in 384 well plates for a modified 5'-RACE RNA-seq library preparation. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Loading the data\n",
    "The first thing we need to do is load the data. In this case, it is provided as an .hdf5 file which, when loaded, represents an AnnData object. This is the format used by the Scanpy library which we'll be using for the majority of this analysis. This object can be thought of as a data matrix in which the rows and columns can be annotated with user-supplied or -generated metadata. A graphical representation is displayed below:\n",
    "\n",
    "![AnnData Description](https://falexwolf.de/img/scanpy/anndata.svg)\n",
    "\n",
    "This format can be used to organize single cell RNA-seq data. The data matrix (.X field) represents the count data itself, where the number of counts for a specific gene are recorded in each column for every cell, represented by a row:\n",
    "\n",
    "| Gene | Cell_1 | Cell_2 | Cell_3 |\n",
    "| ---- | ------ | ------ | -------|\n",
    "| SNORA74 | 1 | 2 | 3 |   \n",
    "| CCDC144B | 1 | 2 | 3 |\n",
    "| LINC00430 | 0 | 1 | 0 |    \n",
    "\n",
    "\n",
    "Individual cell annotions are found in the \".obs\" field. We can add annotations like, what sample a cell belongs to, or what cell type it is, etc. Gene annotations are likewise found in the \".var\" field. Unstructured annotations, those not pertaining to a specific gene or cell but to the whole dataset, are provided in the \".uns\" field. Anndata objects can also contain multiple layers, which are alternatives to the \".X\" field. For instance, when we log normalize our data as part of preprocessing, we can save the original raw count data as a layer for easy access later. The alternative layers can be accessed through the \".layer\" dictionary.\n",
    "\n",
    "Each of these fields, including the \".X\" data matrix, are accessible in the same was as a pandas DataFrames with `.loc` and `.iloc` indexing. Subsetting the AnnData object will also subset its elements accordingly, so an operations like `anndata[list_of_genes,:]` will subset \".X\", \".obs\", and any layers you've created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:02:10.066763Z",
     "start_time": "2020-02-24T19:02:09.493151Z"
    }
   },
   "outputs": [],
   "source": [
    "# Comment and point read_h5ad function towards downloaded dengue dataset.\n",
    "scdata_file = wget.download( \"https://storage.googleapis.com/andersen-lab_temp/dengue_all.hdf5\" )\n",
    "\n",
    "scdata = sc.read_h5ad( scdata_file )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preprocessing\n",
    "\n",
    "The first step in single cell RNA-seq analysis is to prepare them for downstream analysis. We will perform the following steps:\n",
    "\n",
    "1. **Filtration of cells and genes based on quality metrics**   \n",
    "Low quality cells, resulting from damage during processing or incomplete capture, interfere with downstream analyses and must be removed. Genes which provide little information are also removed to speed up computation.  \n",
    "2. **Data normalization**  \n",
    "Normalization eliminates cell/sample-specific biases in capture efficiency and enables comparisons across cells.   \n",
    "3. **Selection of features**  \n",
    "This step selects genes which are variably expressed across cells. By limiting analysis to measurements of these genes we reduce computational requirements and noise from uninteresting or unexpressed genes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quality Control\n",
    "For single cell RNA-seq data, quality control needs to be performed on both cells and genes. To assess the quality of cells we'll use quality measures like the number of molecule counts, the number of genes expressed, and the fraction of counts of mitochondrial genes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:02:11.343092Z",
     "start_time": "2020-02-24T19:02:10.068733Z"
    }
   },
   "outputs": [],
   "source": [
    "# Quickly remove some meaningless annotation\n",
    "scdata = scdata[:,:-6]\n",
    "\n",
    "# Identify mitochondrial genes and spike-in measurements and add their annotation to .var field\n",
    "scdata.var[\"mito\"] = scdata.var_names.str.startswith( \"MT-\" )\n",
    "scdata.var[\"ERCC\"] = scdata.var_names.str.startswith( \"ERCC-\" )\n",
    "\n",
    "# Calculate metrics in regards to mitochondrial reads\n",
    "sc.pp.calculate_qc_metrics( scdata, qc_vars=[ 'mito', 'ERCC' ], inplace=True )\n",
    "\n",
    "print( \"Columns added:\" )\n",
    "for i in scdata.obs_keys():\n",
    "    print( i )\n",
    "\n",
    "print( \"scdata contains {} cells with counts for {} genes\".format( scdata.n_obs, scdata.n_vars ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above code block we calculated all the quality measures we need to perform quality control on our cells. We can see that all together our datasets includes 12263 observations. While the majority of these are likely cells, this number is inflated by the presense of artifacts like dead cells, doublets, and barcode crosstalk. \n",
    "\n",
    "To identify these artifacts, the first thing we'll look at is mitochondrial gene expression. Expression of mitocondrial genes is associated with poor quality cells which have lost cytoplasmic RNA for one reason or another. We'll examine a plot of the fraction of mitochondrial reads vs. number of molecules to determine what the thesholds should be for these quality metrics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:02:11.721137Z",
     "start_time": "2020-02-24T19:02:11.345080Z"
    }
   },
   "outputs": [],
   "source": [
    "# Proposed mitochondrial limit based on graph.\n",
    "proposed_mt_limit = 10\n",
    "\n",
    "# Plot number of fraction of mitochondrial reads vs. molecule counts per cell.\n",
    "ax = plt.subplot()\n",
    "sc.pl.scatter( scdata, \"total_counts\", \"pct_counts_mito\", show=False, ax=ax )\n",
    "plt.axhline( proposed_mt_limit, 0, 60000, linestyle=\"dashed\" )\n",
    "plt.show()\n",
    "\n",
    "# Calculate percentage of cells which would be lost by preposed limit.\n",
    "pct_removed = (scdata[scdata.obs[\"pct_counts_mito\"] > proposed_mt_limit].n_obs / scdata.n_obs ) * 100\n",
    "\n",
    "print( \"{} observations ({:.3f}%) will be removed with proposed mitochondrial limit\".format( scdata[scdata.obs[\"pct_counts_mito\"] > proposed_mt_limit].n_obs, pct_removed ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-04T23:43:40.032002Z",
     "start_time": "2020-02-04T23:43:40.023686Z"
    }
   },
   "source": [
    "In the plot we can see that there is a populations of cells with low molecule counts and a high fraction of mitochondrial gene expression. An appears as if an upper threshold of **10%** mitochondrial gene expression is appropriate so as to remove as much of this population as possible. \n",
    "\n",
    "Next we will look at cell size. We want to filter cells which contain either too little or too many molecule counts or total counts as these would be indicative of dead cells and multiplets, respectively. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:02:13.036055Z",
     "start_time": "2020-02-24T19:02:11.723302Z"
    }
   },
   "outputs": [],
   "source": [
    "# Proposed upper and lower limits based on graphs.\n",
    "lower_counts = 5000\n",
    "upper_counts = 5500000\n",
    "\n",
    "# Plot histogram of molecule counts per cell.\n",
    "sns.distplot( scdata.obs[\"total_counts\"], kde=False )\n",
    "plt.show()\n",
    "\n",
    "# Zoom in on left side of distribution.\n",
    "sns.distplot( scdata.obs.loc[scdata.obs[\"total_counts\"]<100000, \"total_counts\"], kde=False, bins=60 )\n",
    "plt.axvline( lower_counts, 0, 100, linestyle=\"dashed\", color=\"black\" )\n",
    "plt.show()\n",
    "\n",
    "# Zoom in on right side of distribution.\n",
    "sns.distplot( scdata.obs.loc[scdata.obs[\"total_counts\"]>2000000, \"total_counts\"], kde=False, bins=60 )\n",
    "plt.axvline( upper_counts, 0, 20, linestyle=\"dashed\", color=\"black\" )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The histogram of the number of molecule counts per cell is normally distibuted with a peak around 60,000 reads per cell. The upper and lower limits of cell size are selected where the distribution lacks normality.\n",
    "\n",
    "The histrogram of the number of counts per cells reveal a peak of cells with very few counts. Given the distribution of counts, these are likely to be uninformative and likely arise as a a result of debris or dead cells in wells.\n",
    "\n",
    "Likewise, the histogram shows a number of outlier observations which contain massive numbers of counts (>5,500,000 reads). We will remove these as well.\n",
    "\n",
    "We can performed the same analysis on unique gene counts as they are supposed to correlate with molecule counts. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:02:14.861547Z",
     "start_time": "2020-02-24T19:02:13.045259Z"
    }
   },
   "outputs": [],
   "source": [
    "# Proposed lower threshold\n",
    "lower_genes = 650\n",
    "\n",
    "# Calculate the number of unqiue genes expressed by each cell.\n",
    "scdata.obs[\"n_genes_by_counts\"] = [row.count_nonzero() for row in scdata.X]\n",
    "\n",
    "# Plot histogram of unique gene counts.\n",
    "sns.distplot( scdata.obs[\"n_genes_by_counts\"], kde=False )\n",
    "plt.show()\n",
    "\n",
    "# Zoom in on the left side of the distribution.\n",
    "sns.distplot( scdata.obs.loc[scdata.obs[\"n_genes_by_counts\"] < 2000, \"n_genes_by_counts\"], kde=False, bins=60 )\n",
    "plt.axvline( lower_genes, 0, 100, linestyle=\"dashed\", color=\"black\" )\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The histogram of unique gene counts per cell shows that there is a large peak of cells expression less than 500 genes. This peak is made up of dead cells or wells that had only ambient RNA as input. \n",
    "\n",
    "We can now filter the cells based on these limits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:02:16.095251Z",
     "start_time": "2020-02-24T19:02:14.863861Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print( \"Total number of cells: {}\".format( scdata.n_obs ) )\n",
    "\n",
    "scdata = scdata[scdata.obs[\"pct_counts_mito\"] < proposed_mt_limit]\n",
    "print( \"Number of observations after MT filter: {}\".format( scdata.n_obs ) )\n",
    "\n",
    "sc.pp.filter_cells( scdata, min_genes=lower_genes )\n",
    "print( \"Number of observations after min genes filter: {}\".format( scdata.n_obs ) )\n",
    "\n",
    "sc.pp.filter_cells( scdata, min_counts=lower_counts )\n",
    "print( \"Number of observations after min count filter: {}\".format( scdata.n_obs ) )\n",
    "\n",
    "sc.pp.filter_cells( scdata, max_counts=upper_counts )\n",
    "print( \"Number of observations after max count filter: {}\".format( scdata.n_obs ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly, we will filter genes. We only want to included in downstream analyses genes which are expressed in at least 2% of cells. We'll also remove the counts for the ERCC spike-ins. These provide useful metrics for normalization but we don't want to imply anything biological from their values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:02:16.645864Z",
     "start_time": "2020-02-24T19:02:16.096863Z"
    }
   },
   "outputs": [],
   "source": [
    "proposed_gene_limit = 0.02\n",
    "\n",
    "print( \"Total number of genes: {}\".format( scdata.n_vars ) )\n",
    "\n",
    "# Filters out genes.\n",
    "sc.pp.filter_genes( scdata, min_cells=scdata.n_obs * proposed_gene_limit )\n",
    "print( \"Number of genes after cell filter: {}\".format( scdata.n_vars ) )\n",
    "\n",
    "# Remove ERCC spike-ins\n",
    "scdata = scdata[:,~scdata.var[\"ERCC\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normalization\n",
    "Currently, the data we've been working with is represented as counts which are representative of the molecules which were captured in the single cell RNA-seq experiment. Given that the RNA-seq experiment is a random sampling of mRNA in a cell, the number of counts detected in a cell is a function of both the actual number of transcripts in a cell and the capture efficiency of those transcripts. We must then estimate the actual number of transcripts in a cell to be able to compare expression values between cells.  \n",
    "\n",
    "We call the estimated number of transcripts in a cell the *size factor* of a cell. Normalized expression values for a cell can be determine by dividing the measured counts by the size factor of a cell. To calculate these size factors we'll use an approach implemented in *computeSizeFactors* function of the *scran* R package. In brief, this method computes a size factor for clusters of similar cells and then determines cell size factors by deconvoluting the cluster-based size factors. Many alternative methods exist, but a comparison performed by [Büttner et al. 2019](https://doi.org/10.1038/s41592-018-0254-1) found it to perform best.\n",
    "\n",
    "To utilize this method, we'll quickly log normalize a copy of the data and quickly cluster to identify broad cell types. We'll then use these clusters as input into the *computeSizeFactors* function. Later on in the tutorial we'll go into these steps in more depth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:02:30.647722Z",
     "start_time": "2020-02-24T19:02:16.649133Z"
    }
   },
   "outputs": [],
   "source": [
    "# Create a copy of the data\n",
    "scdata_copy = scdata.copy()\n",
    "\n",
    "# Log normalize data\n",
    "sc.pp.normalize_per_cell( scdata_copy, counts_per_cell_after=1e6, key_n_counts=\"total_counts\" )\n",
    "sc.pp.log1p( scdata_copy )\n",
    "\n",
    "# Quickly cluster the data\n",
    "sc.pp.pca( scdata_copy, svd_solver=\"arpack\", n_comps=15 )\n",
    "sc.neighbors.neighbors( scdata_copy, n_pcs=15, use_rep=\"X_pca\", random_state=42 )\n",
    "sc.tl.louvain( scdata_copy, key_added=\"groups\", resolution=0.5 )\n",
    "\n",
    "# Collect cluster assignments and data matrix for input into R.\n",
    "input_groups = scdata_copy.obs[\"groups\"]\n",
    "data_mat = scdata.X.T.toarray()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because *scran* is an R package, we need to use Rpy2 to interface with R. This block takes the cluster asignments and the transposed data matrix as input and returns a list of size factors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:04:45.935558Z",
     "start_time": "2020-02-24T19:02:30.649927Z"
    }
   },
   "outputs": [],
   "source": [
    "%%R -i data_mat -i input_groups -o size_factors\n",
    "library( scran )\n",
    "\n",
    "size_factors = computeSumFactors( data_mat, clusters=input_groups, min.mean=0.1, positive=TRUE )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:04:47.147295Z",
     "start_time": "2020-02-24T19:04:45.937087Z"
    }
   },
   "outputs": [],
   "source": [
    "# Delete the data copy to free up space\n",
    "del scdata_copy\n",
    "\n",
    "# Visualize the estimated size factors\n",
    "scdata.obs[\"size_factors\"] = size_factors\n",
    "\n",
    "sc.pl.scatter( scdata, \"size_factors\", \"total_counts\" )\n",
    "sc.pl.scatter( scdata, \"size_factors\", \"n_genes_by_counts\" )\n",
    "\n",
    "sns.distplot( size_factors, bins=50, kde=False )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By plotting the size factors we can see that they mostly correlate with the molecule counts and number of unique genes.\n",
    "\n",
    "Next we'll save the raw count data in a new layer for access later, normalize the data matrix by dividing by the size factors, and log transform the data. This transformation is necessary to identify strong relative differences between gene expression. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:04:52.121250Z",
     "start_time": "2020-02-24T19:04:47.149120Z"
    }
   },
   "outputs": [],
   "source": [
    "# Keep the count data in a counts layer\n",
    "scdata.layers[\"counts\"] = scdata.X.copy()\n",
    "\n",
    "# Normalize and log tranform data matrix.\n",
    "scdata.X /= scdata.obs[\"size_factors\"].values[:,None]\n",
    "sc.pp.log1p( scdata )\n",
    "\n",
    "# Convert matrix to sparse matrix to free up space. Sparse matrix doesn't store zero counts.\n",
    "scdata.X = sp.sparse.csr_matrix( scdata.X )\n",
    "\n",
    "# Store the full data set in 'raw' as log-normalised data for statistical testing\n",
    "scdata.raw = scdata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Batch Correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As discussed earlier, the dataset contains cells from six different individuals. These samples were processed and sequenced individually so we need to identify and remove batch effects from the data. Batch effects and biological effects are often overlapping so any type of batch effect correction will necessarily remove some biological structure from the data. For this reason, we use the batch corrected data for clustering, but use the non-batch corrected data for calculating differential expression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:05:28.907969Z",
     "start_time": "2020-02-24T19:04:52.122964Z"
    }
   },
   "outputs": [],
   "source": [
    "sc.pp.combat( scdata, key=\"PatientID\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Highly Variable Genes\n",
    "The clustering and dimensionality reduction which we will perform compares cells based on their gene expression profiles. However, not all genes provide useful information towards this end. Some are lowly expressed and their counts when present amount to noise in the data, while other provide interesting biological signal. We will remove the former group of genes and retain the latter to not only reduce the size of the dataset, but to highlight relevent biological structure.\n",
    "\n",
    "To identify the useful set of genes we select from the genes with the highest variability; genes which are highly expressed in some cells, and lowly expressed in others. Here we use a technique of identifying highly variable genes from the Chromium 10x preprocessing software *CellRanger*. The number of genes to retain is an important decision. Enough genes need to be retains so as not to lose any biological signal, but not too many genes that noise is included from irrelevent genes. A number of strategies can be employed to determine the optimal amount (see [here](https://osca.bioconductor.org/feature-selection.html#hvg-selection)), but for our purposes we'll just select the top *n* genes. Typically between 1000 and 5000 genes are selected and for our analysis we'll include 4000. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:05:37.824747Z",
     "start_time": "2020-02-24T19:05:34.251620Z"
    }
   },
   "outputs": [],
   "source": [
    "sc.pp.highly_variable_genes( scdata, flavor='cell_ranger', n_top_genes=4000 )\n",
    "print( \"Number of highly variable genes: {}\".format( sum( scdata.var[\"highly_variable\"] ) ) )\n",
    "\n",
    "sc.pl.highly_variable_genes( scdata )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dimensionality Reduction\n",
    "When analyzing single cell RNA-seq data, we aim to compare cells based on their expression across genes, i.e. to identify distinct transcriptomic profiles. However, expression of different genes is often correlated if they are affected by the same underlying biological process. In order to study these underlying biological processes, its necessary to compress the information from individual genes into a single dimension. Dimensionality reduction aims to reduce the number of dimensions needed to represent a dataset and presents a number of advantages:\n",
    "* Reduces the computation load of downstream analyses.\n",
    "* Reduced noise in dataset by averaging multiple correlated genes.\n",
    "* Enables visualization of high dimensional dataset.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Principal Component Analysis\n",
    "The first technique we'll be using is Principal Component Analysis, or PCA. PCA reduces dimensionality by finding axes in the high dimensional data which capture the largest amount of variation. \n",
    "\n",
    "To perform PCA, an axis through the high dimensional data, also called a principal component (PC), is found which captures the greatest variance across cells. Next, another component is selected that is perpendicular to the first and captures the greatest remaining variance. This process can be repeated until a set number of principal components are found. \n",
    "\n",
    "Here we'll perform PCA on the normalized data and compute the first 50 principal components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:05:51.727579Z",
     "start_time": "2020-02-24T19:05:37.879338Z"
    }
   },
   "outputs": [],
   "source": [
    "# Calculate the visualizations\n",
    "sc.pp.pca( scdata, n_comps=50, use_highly_variable=True, random_state=42 )\n",
    "sc.pl.pca( scdata, color=\"n_counts\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this plot we can already see distinct populations of cells, but, given that we've only plotted two principal components, its not that informative beyond that. One thing we can look at is the genes which contribute significantly to the first two principal components and ask whether they represent anything meaningful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:05:53.077573Z",
     "start_time": "2020-02-24T19:05:51.729537Z"
    }
   },
   "outputs": [],
   "source": [
    "sc.pl.pca_loadings( scdata, components=[1,2], include_lowest=False )\n",
    "sc.pl.pca( scdata, color=[\"S100A9\", \"NKG7\"] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We computed 50 principal components which represents a rather compact version of the dataset, though this is an arbitrary number. One question we can ask is how many principal components are neccessary to include for downstream analyses. Including more principal components means including more biological signal, but at the cost of more noise.\n",
    "\n",
    "One heuristic we can use to select an appropriate number of principal components is an elbow plot, a plot which ranks principical components based on the percentage of variance they explain. We then include all principal components which lie above an elbow point on the plot where components above the point explain much more variance than the remaining components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:05:53.490725Z",
     "start_time": "2020-02-24T19:05:53.079493Z"
    }
   },
   "outputs": [],
   "source": [
    "sc.pl.pca_variance_ratio( scdata, log=True, n_pcs=50 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By looking at the graph we can see that we should probably include the top 17 principal components. Feel free to play around with the number of principal components to include and the effect it has on later analyses. A general rule is that true biological structure will be independent of the number of principal components included in clustering."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Clustering\n",
    "We'll now algorithmically define the populations of cells that we saw in the PCA plot with unsupervised clustering. \n",
    "\n",
    "The type of clustering we'll be performing is called graph-based clustering. First, we build a neighborhood graph where each node is a cell with edges connecting it to its nearest neighbors in terms of transcriptional profiles. Edges are then weighted based on the similarily to the neighboring cells. This graph is then used as input to a clustering algorithm which identifies clusters as cells which are more connected to eachother than to cells of other clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:06:46.527317Z",
     "start_time": "2020-02-24T19:06:44.753898Z"
    }
   },
   "outputs": [],
   "source": [
    "sc.pp.neighbors( scdata, n_pcs=17, random_state=42 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We next need to decide which clustering algorithm to use. There are essentially two commonly used algorithms, Louvain and Leiden. Both iteratively group cells together and attempt to maximize the difference between the number of edges in a cluster and the number of edges expected to arise due to chance (modularity). We'll perform both and compare results to determine which is better for this dataset.\n",
    "\n",
    "One other factor we need to consider is the resolution parameter which both algorithms require. This represents the granularity of the resulting clusters. Higher values generate a greater number of clusters and typically values of between 0.4 and 1.5 are seen depending on the size of the dataset. For our purposes, and since we only need to cluster broad cell types now, a resolution of 0.5 will be used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:06:54.625490Z",
     "start_time": "2020-02-24T19:06:51.897332Z"
    }
   },
   "outputs": [],
   "source": [
    "sc.tl.leiden( scdata, resolution=0.5, key_added=\"leiden_r0.5\", random_state=42 )\n",
    "sc.tl.louvain( scdata, resolution=0.5, key_added=\"louvain_r0.5\", random_state=42 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To visualize the clustering results, we apply further dimensional reduction using a non-linear, unlike PCA, technique called UMAP. Similar to t-SNE, the UMAP algorithm attempts to accurately convey the local and global distances between cells' transcriptomic profiles in two dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:07:10.435505Z",
     "start_time": "2020-02-24T19:06:56.780273Z"
    }
   },
   "outputs": [],
   "source": [
    "sc.tl.umap( scdata, random_state=42 )\n",
    "sc.pl.umap( scdata, color=[\"louvain_r0.5\", \"leiden_r0.5\"] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The clustering results are fairly comparable and that is probably a result of the clustering algorithms being fairly similar. We'll use the leiden algorithm for the remainder of this tutorial as there are some arguments in the literature that it maintains well-connected communities better than louvain clustering (see [Traag et al. 2018](https://arxiv.org/ct?url=https%3A%2F%2Fdx.doi.org%2F10.1038%2Fs41598-019-41695-z&v=1b4d93af) for more information).\n",
    "\n",
    "Feel free to compare the results with the louvain algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# T Cell Selection\n",
    "To determine which of the clusters we just calculated correspond to T cells, we're going to look at the expression of canonical T cell markers, CD3 and the T cell receptor alpha constant chain, in these clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:07:24.175654Z",
     "start_time": "2020-02-24T19:07:23.624850Z"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure( dpi=100, figsize=(10,10) )\n",
    "ax = plt.subplot()\n",
    "sc.pl.umap( scdata, color=\"leiden_r0.5\", ax=ax, legend_loc=\"on data\" )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:07:38.522331Z",
     "start_time": "2020-02-24T19:07:37.462180Z"
    }
   },
   "outputs": [],
   "source": [
    "sc.pl.violin( scdata, [\"CD3E\", \"TRAC\"], groupby=\"leiden_r0.5\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note: This may not be deterministic. Continue on your own, or have a look at the \"Workshop_Notebook_complete\" notebook for my results, if yours look different.**\n",
    "\n",
    "By using the UMAP embeddings and the expression of canonical T cell markers, we can observe that T cells are likely found in clusters 3, 4, and 13 \n",
    "\n",
    "Cluster 8 and 13 are also a candidate for containing T cells, but we won't include them in the T cell group. Both express some level of CD3 and TRAC, but that likely is a result of multiplets or other technical phenomenons. The other puntative T cell clusters group closely together in the UMAP embedding while putative T cells in clusters 8 and 13 are intermediate to larger non-T cell clusters and T cells.\n",
    "\n",
    "It is actually likely that all clusters contain some percentage of multiplet cells and cells which are incorrectly clusters. The next step in our analysis will be to identify cells in our putative T cell clusters that fit either of these cases using a method developed by [Szabo et al. 2019](https://doi.org/10.1038/s41467-019-12464-3). In brief, we generate a blacklist of genes that are overexpressed in non-T cell cluster and use those to remove contaminating cells in the T cell cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:08:03.698211Z",
     "start_time": "2020-02-24T19:08:03.687413Z"
    }
   },
   "outputs": [],
   "source": [
    "# Merge and rename T cell cluster\n",
    "scdata.obs[\"cluster\"] = (scdata.obs[\"leiden_r0.5\"]\n",
    "                         .map(lambda x: {\"3\" : \"T Cell\", \"4\" : \"T Cell\", \"13\" : \"T Cell\"}.get(x, x) )\n",
    "                         .astype(\"category\") )\n",
    "scdata.obs[\"cluster\"].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate Blacklist\n",
    "To generate the blacklist of genes, we conduct differential expression analysis pairwise between the T cell cluster and each non-T cell cluster. We then add all genes have an adjusted p-value of greater than 0.001 and greater than 5-fold enrichment in non-T cell clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:09:18.493753Z",
     "start_time": "2020-02-24T19:08:27.751590Z"
    }
   },
   "outputs": [],
   "source": [
    "# Perform differential expression analysis, comparing each cluster with the putative T cell cluster.\n",
    "sc.tl.rank_genes_groups( scdata, groupby=\"cluster\", groups=[0,1,2,5,6,7,8,9,10,11,12,14,15], reference=\"T Cell\", method=\"wilcoxon\", use_raw=False, n_genes=-1 )\n",
    "\n",
    "# Add differentially expressed genes for each cluster to the blacklist.\n",
    "blacklist = list()\n",
    "for i in [\"0\",\"1\",\"2\",\"5\",\"6\",\"7\",\"8\",\"9\",\"10\",\"11\",\"12\",\"14\",\"15\"]:\n",
    "    blacklist_genes = (scdata.uns[\"rank_genes_groups\"][\"logfoldchanges\"][i] > 5 ) & ( scdata.uns[\"rank_genes_groups\"][\"pvals_adj\"][i] < 0.001 )\n",
    "    blacklist.extend( scdata.uns[\"rank_genes_groups\"][\"names\"][i][blacklist_genes] )\n",
    "blacklist = np.unique( blacklist )\n",
    "print( \"{} genes added to the blacklist\".format( len( blacklist ) ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate Whitelist \n",
    "We also generate a whitelist of genes which are positively enriched in any T cell cluster. In this way we avoid inclusion of genes which are differentially expressed in a limited number of T cells. \n",
    "\n",
    "We perform differential expression analysis again, but this time comparing all putative T cell clusters to all other identified clusters. We use the same requirements as noted for generating the gene blacklist."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:09:56.685399Z",
     "start_time": "2020-02-24T19:09:40.194198Z"
    }
   },
   "outputs": [],
   "source": [
    "# Perform differential expression analysis, comparing each of the T cell clusters to the rest of the clusters.\n",
    "sc.tl.rank_genes_groups( scdata, groupby=\"leiden_r0.5\", groups=[3,4,13], reference=\"rest\", method=\"wilcoxon\", use_raw=False, n_genes=-1 )\n",
    "\n",
    "# Combine differentially expressed genes\n",
    "whitelist = list()\n",
    "for i in [\"3\", \"4\", \"13\"]:\n",
    "    whitelist_genes = (scdata.uns[\"rank_genes_groups\"][\"logfoldchanges\"][i] > 5 ) & ( scdata.uns[\"rank_genes_groups\"][\"pvals_adj\"][i] < 0.001 )\n",
    "    whitelist.extend( scdata.uns[\"rank_genes_groups\"][\"names\"][i][whitelist_genes] )\n",
    "whitelist = np.unique( whitelist )\n",
    "\n",
    "print( \"{} genes added to the whitelist\".format( len( whitelist ) ) )\n",
    "print( whitelist )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:11:22.449486Z",
     "start_time": "2020-02-24T19:11:22.439770Z"
    }
   },
   "outputs": [],
   "source": [
    "blacklist = [x for x in blacklist if x not in whitelist]\n",
    "print( \"{} genes ultimately in blacklist\".format( len( blacklist ) ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we would expect, the T cell whitelist genes include those which we would typically call T cell marker genes. Broad T cell markers like *CD3*, *LTB*, and T cell receptor genes, cytotoxic markers like *GZMA* and *GZMK*, and cytokines like *CCL5*, *CCR7*, and *IL32*.\n",
    "\n",
    "Thankfully, the whitelist only shared 4 genes with the blacklist."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remove Non-T cells\n",
    "We'll now use the blacklist to remove cells which have been erroneously included in T cell clusters. We will then remove cells which express the blacklist genes above a certain threshold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:11:51.117934Z",
     "start_time": "2020-02-24T19:11:46.346155Z"
    }
   },
   "outputs": [],
   "source": [
    "# Calculate the percentage of normalized counts for blacklisted genes detected for each cell.\n",
    "scdata.var[\"blacklist\"] = scdata.var_names.isin( blacklist )\n",
    "temp = sc.pp.calculate_qc_metrics( scdata, qc_vars=[\"blacklist\"] )\n",
    "scdata.obs[\"pct_blacklist\"] = temp[0][\"pct_counts_blacklist\"]\n",
    "\n",
    "# Calculate mean and absolute deviation of expression of blacklisted genes in T cell cluster.\n",
    "mean_pct_blacklist = scdata.obs.loc[(scdata.obs[\"cluster\"]==\"T Cell\"), \"pct_blacklist\"].mean()\n",
    "mad_pct_blacklist = scdata.obs.loc[(scdata.obs[\"cluster\"]==\"T Cell\"), \"pct_blacklist\"].mad()\n",
    "\n",
    "# Plot distribution of blacklisted gene expression\n",
    "plt.figure( dpi=100 )\n",
    "ax = plt.subplot()\n",
    "sns.distplot( scdata.obs.loc[(scdata.obs[\"cluster\"]==\"T Cell\"), \"pct_blacklist\"], bins=np.arange(1,20,0.5), label=\"Cells in T cell Clusters\", kde=False )\n",
    "sns.distplot( scdata.obs.loc[(scdata.obs[\"cluster\"]!=\"T Cell\"), \"pct_blacklist\"], bins=np.arange(1,20,0.5), label=\"Cells in Non-T cell Clusters\", kde=False )\n",
    "plt.axvline( mean_pct_blacklist + 2.5*mad_pct_blacklist, 0, 20, linestyle=\"dashed\", color=\"black\" )\n",
    "plt.ylabel( \"Cells\" )\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "blacklist_threshold = mean_pct_blacklist + 2.5 * mad_pct_blacklist\n",
    "print( \"Threshold determined to be: {:.2f}%\".format( blacklist_threshold ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interestingly, and likely because of the high per-cell sequencing depth of this dataset, we see a certain amount of expression of the blacklisted genes in every cell. However, we can exploit the bi/tri-modal nature of the distribution to adequately filter T cells. Szabo et al. ascribes the expression of blacklisted genes in the T cell cluster to PCR recombination, where a molecule originating from one cell receives the cell-identifying barcode of a different cell at an intermediate point in PCR. \n",
    "\n",
    "We can see that the T cell distribution appears to be normally distributed with a slight positive skew. A simple threshold would then be anywhere from 2-3 mean absolute deviations above the mean of the T cell distribution. We use the mean absolute deviation here instead of the standard deviation because its less sensitive to outliers which this dataset doesn't have but others typically do. We consider any cell with expression of the blacklisted genes above this threshold and any cell in a non-T cell cluster to be a non-T cell and remove them from further analysis (Plot the dataset and color by cells below this threshold to see why we don't include all cells with expression of blacklisted genes below this threshold)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:12:32.826654Z",
     "start_time": "2020-02-24T19:12:32.760082Z"
    }
   },
   "outputs": [],
   "source": [
    "T_before = scdata.obs[\"cluster\"].value_counts()[\"T Cell\"]\n",
    "scdata = scdata[(scdata.obs[\"cluster\"]==\"T Cell\") & (scdata.obs[\"pct_blacklist\"]<blacklist_threshold)]\n",
    "T_after = scdata.n_obs\n",
    "\n",
    "print( \"Number of cells removed by blacklist gene filter: {}\".format( T_before - T_after ) )\n",
    "print( \"Number of T cells identified: {}\".format( T_after ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# T Cell Clustering\n",
    "Now that T cells have been identified, we'll cluster again to determine what populations we can identify. We will use the same procedure as before, calculating highly variable genes, dimensionality reduction, followed by clustering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:12:36.605264Z",
     "start_time": "2020-02-24T19:12:32.828671Z"
    }
   },
   "outputs": [],
   "source": [
    "# Identify highly variable genes\n",
    "sc.pp.highly_variable_genes( scdata, flavor='cell_ranger', n_top_genes=2000 )\n",
    "print( \"Number of highly variable genes: {}\".format( sum( scdata.var[\"highly_variable\"] ) ) )\n",
    "\n",
    "# Dimensionality reduction using PCA\n",
    "sc.pp.pca( scdata, n_comps=50, use_highly_variable=True, random_state=42 )\n",
    "sc.pl.pca_variance_ratio( scdata, log=True, n_pcs=50 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We lose a bit of structure, biological or otherwise, when limiting our analysis to T cells, so we only need to include the top 6 principal components for clustering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:12:39.774714Z",
     "start_time": "2020-02-24T19:12:36.607511Z"
    }
   },
   "outputs": [],
   "source": [
    "# Build neighborhood graph\n",
    "sc.pp.neighbors( scdata, n_pcs=6, random_state=42 )\n",
    "\n",
    "# Cluster at a middling resolution\n",
    "sc.tl.leiden( scdata, resolution=0.4, key_added=\"subset\", random_state=42 )\n",
    "\n",
    "# Calculate UMAP embedding\n",
    "sc.tl.umap( scdata, random_state=42 )\n",
    "\n",
    "# \n",
    "plt.figure( dpi=100, figsize=(6,6) )\n",
    "ax = plt.subplot()\n",
    "sc.pl.umap( scdata, color=[\"subset\"], ax=ax, title=\"T Cell Clustering: Resolution 0.4\" )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clustering has identified five T cell clusters. Playing around with the resolution parameter of the leiden function may find more or less clusters, though these five appear to be the most stable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# T Cell Subset Characterization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CD4 vs. CD8\n",
    "The ultimate differentiator of T cells is whether they're CD4- or CD8-positive. We can visually assess this using the UMAP embedding and color each cell by the ratio of its expression of CD4 to CD8. Because drop-out, where a transcript is randomly not sampled during the RNA-seq experiment, is quite common we'll color cells by the mean ratio of CD4 and CD8 in the cluster it has been assigned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:12:40.190055Z",
     "start_time": "2020-02-24T19:12:39.776713Z"
    }
   },
   "outputs": [],
   "source": [
    "# Calculate the log ratio of CD4 to CD8 expression\n",
    "scdata.obs[\"CD4\"] = scdata[:,\"CD4\"].X\n",
    "scdata.obs[\"CD8\"] = scdata[:,\"CD8A\"].X\n",
    "scdata.obs[\"CD4/CD8\"] = scdata.obs.groupby(\"subset\")[\"CD4\"].transform( \"mean\" ) - scdata.obs.groupby( \"subset\" )[\"CD8\"].transform( \"mean\" )\n",
    "scdata.obs = scdata.obs.drop( columns=[\"CD4\", \"CD8\"] )\n",
    "\n",
    "plt.figure( dpi=100, figsize=(6,6) )\n",
    "ax = plt.subplot()\n",
    "sc.pl.umap( scdata, color=[\"CD4/CD8\"], ax=ax, vmin=-3.5, vmax=3.5, color_map=\"seismic\" )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clusters 0 and 3 are likely CD4-positive T cells, while clusters 1, 2, and 4 are CD8-positive T cells."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Subset Identification\n",
    "To give sense to the identified subsets, besides merely CD4 and CD8 expression, we need to determine the genes whose expressions differentiates the clusters. We can then use these \"marker\" genes to assign clusters to canonical cell types and determine functional differences between clusters.\n",
    "\n",
    "Marker genes are identified by their differential expression between clusters. Genes which are more strongly differentially expression are more likely to cause seperation between clusters of cells. To quantify differential expression we can use any number of statistical tests (t-test, wilcoxon test, linear regression) to make pairwise comparisons between clusters. For our purposes, we'll use the Wilcoxon rank sum test which is comparable to more complex methods. For a thorough comparison of differential expression methods, check out [Soneson and Robinson 2018](https://doi.org/10.1038/nmeth.4612)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:12:43.813100Z",
     "start_time": "2020-02-24T19:12:40.191652Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Perform DE analysis\n",
    "sc.tl.rank_genes_groups( scdata, groupby=\"subset\", method=\"wilcoxon\", use_raw=True, n_genes=100 )\n",
    "\n",
    "# Plot top 20 differentially expressed genes for each cluster.\n",
    "sc.pl.rank_genes_groups( scdata, n_genes=20, ncols=3, sharey=False )\n",
    "\n",
    "# Gather all differentially expressed genes into a dataframe for easy access.\n",
    "marker_genes = pd.DataFrame( scdata.uns[\"rank_genes_groups\"][\"names\"] )\n",
    "marker_genes.head(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use the differentially expressed genes to determine what biological state is represented by each of the clusters. This isn't entirely straightforward and is likely the most challenging part of scRNA-seq analysis. An optimal identification method would be to look at the overlap of differentially expressed genes with a list of known markers genes, but often these aren't available for the more novel and less defined cell subsets we're looking at here. \n",
    "\n",
    "Below I've laid out the rational for how I would annotate the clusters.\n",
    "\n",
    "Putative cluster labels:\n",
    "- Cluster 0 - **CD4-C1-Resting** - Express markers of central memory and naive phenotype (*CCR7*, *SELL*), and resting (*LEF1*). Also express ribosomal genes indicating active translation.\n",
    "- Cluster 1 - **CD8-C1-NKGY** - Cytotoxic T cell subset. Expresses markers of cytotoxicity (*PRF1*, *NKG7*, *GNLY*, *GZMH*), exhaustion (*CST7*). \n",
    "- Cluster 2 - **CD8-C2-KLRB1** - Cytotoxic T cell subset. Expresses markers of cytotoxicity (*KLRB1*, *KLRG1*, *NKG7*, *GZMK*, *PRF1*)\n",
    "- Cluster 3 - **CD4-C2-Activated** - Activated T cell subset. Express markers of T cell motility, proliferation, and development (*COTL1*, *ADAM19*, *TRADD*, *GSTK1*, *ITGB1*)\n",
    "- Cluster 4 - **CD8-C3-Resting** - Express markers of central memory and naive phenotype, and resting (*CCR7*)\n",
    "\n",
    "We'll annotate the clusters so that downstream analysis is easier to read."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:12:44.101138Z",
     "start_time": "2020-02-24T19:12:43.815608Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Rename clusters based on assignment above.\n",
    "scdata.rename_categories( \"subset\", [\"CD4-C1-Resting\",\"CD8-C1-NKGY\",\"CD8-C2-KLRB1\",\"CD4-C2-Activated\",\"CD8-C3-Resting\"] )\n",
    "scdata.obs[\"subset\"].value_counts()\n",
    "\n",
    "# Plot cluster assignments on UMAP embeddings.\n",
    "plt.figure( dpi=100, figsize=(6,6) )\n",
    "ax=plt.subplot()\n",
    "sc.pl.umap( scdata, color=\"subset\", legend_loc=\"on data\", show=False, ax=ax )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One other way we can look at these marker genes is through a heatmap of the highest differentially expressed genes. This will show us which clusters are most similar and how their markers genes are expressed in other clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:12:45.395457Z",
     "start_time": "2020-02-24T19:12:44.102573Z"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure( dpi=200 )\n",
    "sc.pl.rank_genes_groups_heatmap( scdata, key=\"rank_genes_groups\", n_genes=20, show=False, figsize=(15,15), standard_scale=\"var\" )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to the umap embeddings, T cell clusters mostly group together based on their CD4-/CD8-positive character and their activation state.\n",
    "\n",
    "We can also see that differential gene expression analysis had difficulty determining marker genes for the resting CD4-/CD8-positive populations. Likely, these subsets are defined based on their lack of expression of certain activation markers, which we wouldn't see in this or the previous charts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# T Cell Functional Comparison\n",
    "The dataset we're using allows us to determine T cell-specific effects which differentiate acute and severe dengue fever. The disease status, whether an individual is healthy or has acute or sever dengue fever, is stored for each cell in the .obs field. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compositional Analysis\n",
    "An initial comparison we can perform is a compositional one where we look at whether infection causes a significant difference in the size of various T cell populations. It should be noted that this type of analysis requires a sufficient number of cells to accurately asses populations sizes and account for expected variability. Unfortunately there are no dedicated tools to perform this type of analysis, so visual comparison of the composition is as informative as we can be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:12:45.857249Z",
     "start_time": "2020-02-24T19:12:45.398205Z"
    }
   },
   "outputs": [],
   "source": [
    "# Create a cross table comparing cluster size to disease status. \n",
    "cells_patient = pd.crosstab( scdata.obs[\"Disease_Status\"], scdata.obs[\"subset\"] ).apply( lambda r: np.round( r/r.sum()*100, 2 ), axis=1 )\n",
    "\n",
    "# Plot the cross table as a heatmap.\n",
    "plt.figure( dpi=100, figsize=(7,7) )\n",
    "ax=plt.subplot()\n",
    "sns.heatmap( cells_patient, square=True, annot=True, fmt=\".1f\", cmap=\"Blues\", ax=ax )\n",
    "plt.xlabel( \"Cluster\" )\n",
    "plt.ylabel( \"Disease Status\" )\n",
    "for t in ax.texts: t.set_text( t.get_text() + \"%\" )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above cross table we can see that patients with acute dengue fever have a higher proportion of resting CD4-positive T cells while patients with severe dengue fever have a high proportion of activated CD4-positive T cells. Acute dengue fever patients also have elevetated levels of one of cytotoxic CD8-positive T cell clusters, while individuals with severe dengue fail to see any elevation of these populations above control individual levels. \n",
    "\n",
    "To view this in another way, we can visualize shifts in cellular densities between disease statuses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:12:46.641471Z",
     "start_time": "2020-02-24T19:12:45.859484Z"
    }
   },
   "outputs": [],
   "source": [
    "# Calculate cellular densities for each of the disease statuses.\n",
    "sc.tl.embedding_density( scdata, groupby=\"Disease_Status\", key_added=\"umap_density_disease_status\" )\n",
    "\n",
    "# Plot the cellular density for each of the disease statuses.\n",
    "plt.figure( dpi=200 )\n",
    "sc.pl.embedding_density( scdata, key=\"umap_density_disease_status\", group=\"control\", color_map=\"Blues\" )\n",
    "sc.pl.embedding_density( scdata, key=\"umap_density_disease_status\", group=\"dengue fever\", color_map=\"Blues\" )\n",
    "sc.pl.embedding_density( scdata, key=\"umap_density_disease_status\", group=\"severe dengue\", color_map=\"Blues\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visually, there is a striking difference in the cellular densities. Most clusters are evenly represented in the healthy control individuals, while individuals with acute dengue and severe dengue show much higher proportions of activated subsets. Individuals with acute dengue fever fail to have an activated CD4-postive populations as large individuals with severe dengue fever. Both acute and severe dengue fever are associated with large populations of cytotoxic T cells, but these represented distinct populations. Additionally, because the size of this population in individuals with severe dengue is not visually different from the population in healthy individuals, it would seem that severe dengue correlates with the non-presense of suitable cytotoxic T cell populations for virus clearance.\n",
    "\n",
    "This all seems to align with current literature where it has been suggested that both CD4- and CD8-postive T cell populations are import for clearance of dengue virus and progression to severe dengue results from immopathology as a result of the expansion of pre-existing crossreactive and low affinity memory T cell populations against the virus.\n",
    "\n",
    "With this dataset we can also test the theory that the immunopathology associated with severe dengue results from the excess production of inflammatory cytokines. We'll look for evidence of their production in differentially expressed genes between acute and severe dengue fever. \n",
    "\n",
    "**Note**: Differentially expression between conditions is typically performed within clusters. Thus, we'll only look at differentially expressed genes within clusters which are likely to be producing cytokines, the activated CD4-postive and one of the activated CD8-positive populations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:12:48.114603Z",
     "start_time": "2020-02-24T19:12:46.643108Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Subset dataset\n",
    "cd8_activated = scdata[scdata.obs[\"subset\"]==\"CD8-C2-KLRB1\"]\n",
    "cd4_activated = scdata[scdata.obs[\"subset\"]==\"CD4-C2-Activated\"]\n",
    "\n",
    "# Filter out genes which aren't expressed to reduced the number of hypothesies tested.\n",
    "sc.pp.filter_genes( cd8_activated, min_cells=1 )\n",
    "sc.pp.filter_genes( cd4_activated, min_cells=1 )\n",
    "\n",
    "# Perform differential gene analysis\n",
    "sc.tl.rank_genes_groups( cd8_activated, \n",
    "                         groupby=\"Disease_Status\", \n",
    "                         groups=[\"severe dengue\"], \n",
    "                         reference=\"dengue fever\", \n",
    "                         key_added=\"severe_vs_acute\",\n",
    "                         method=\"wilcoxon\",\n",
    "                         n_genes=-1 )\n",
    "sc.tl.rank_genes_groups( cd4_activated, \n",
    "                         groupby=\"Disease_Status\", \n",
    "                         groups=[\"severe dengue\"], \n",
    "                         reference=\"dengue fever\", \n",
    "                         key_added=\"severe_vs_acute\",\n",
    "                         method=\"wilcoxon\",\n",
    "                         n_genes=-1 )\n",
    "\n",
    "# Construct dataframe from results\n",
    "cd4_results = cd4_activated.uns[\"severe_vs_acute\"]\n",
    "cd8_results = cd8_activated.uns['severe_vs_acute']\n",
    "cd8_df = pd.DataFrame( { \"names\": cd8_results[\"names\"][\"severe dengue\"],\n",
    "                         \"scores\": cd8_results[\"scores\"][\"severe dengue\"],\n",
    "                         \"logfoldchanges\": cd8_results[\"logfoldchanges\"][\"severe dengue\"], \n",
    "                         \"pvals\": cd8_results[\"pvals\"][\"severe dengue\"],\n",
    "                         \"pvals_adj\": cd8_results[\"pvals_adj\"][\"severe dengue\"] } )\n",
    "cd4_df = pd.DataFrame( { \"names\": cd4_results[\"names\"][\"severe dengue\"],\n",
    "                         \"scores\": cd4_results[\"scores\"][\"severe dengue\"],\n",
    "                         \"logfoldchanges\": cd4_results[\"logfoldchanges\"][\"severe dengue\"], \n",
    "                         \"pvals\": cd4_results[\"pvals\"][\"severe dengue\"],\n",
    "                         \"pvals_adj\": cd4_results[\"pvals_adj\"][\"severe dengue\"] } )\n",
    "\n",
    "# Remove outlier genes\n",
    "cd8_df[\"neglogpval_adj\"] = -np.log( cd8_df[\"pvals_adj\"] )\n",
    "cd4_df[\"neglogpval_adj\"] = -np.log( cd4_df[\"pvals_adj\"] )\n",
    "cd8_df = cd8_df.loc[cd8_df[\"logfoldchanges\"]>0]\n",
    "cd4_df = cd4_df.loc[cd4_df[\"logfoldchanges\"]>0]\n",
    "cd8_df = cd8_df.loc[cd8_df[\"logfoldchanges\"]<15]\n",
    "cd4_df = cd4_df.loc[cd4_df[\"logfoldchanges\"]<15]\n",
    "\n",
    "# Pull out genes to label\n",
    "labels = [\"IFI6\", \"IFITM3\", \"IFIT3\", \"IFIT2\", \"IL6ST\", \"ISG15\", \"ISG20\", \"STAT2\", \"STAT1\", \"OAS3\", \"OAS1\" ]\n",
    "labels_df = cd8_df.loc[cd8_df[\"names\"].isin( labels )]\n",
    "labels_df[\"logfoldchanges\"] = labels_df[\"logfoldchanges\"] + 0.1\n",
    "labels_df[\"neglogpval_adj\"] = labels_df[\"neglogpval_adj\"] + 1.25\n",
    "\n",
    "# Perform the plotting\n",
    "plt.figure( dpi=100, figsize=(5,5) )\n",
    "ax=plt.subplot()\n",
    "sns.scatterplot( cd8_df[\"logfoldchanges\"], cd8_df[\"neglogpval_adj\"], ax=ax )\n",
    "plt.xlabel( \"logFC\")\n",
    "plt.ylabel( \"-logP-value\")\n",
    "plt.title( \"CD8-C2-KLRB1: Severe vs. Acute Dengue\")\n",
    "for i in labels_df.iterrows():\n",
    "    ax.annotate( i[1][\"names\"], (i[1][\"logfoldchanges\"], i[1][\"neglogpval_adj\"]) )\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the activated CD8 cells, we don't see increased expression of any cytokines besides IL6. This isn't entirely unexpected because cytokine transcripts are produced in very few copies even when upregulated, and might not be present in the circulating T cell populations we've sampled. Instead, we can look at indirect measures of cytokine production. Genes related to the interferon response, including interferon inducible proteins genes and interferon stimulated genes, are significantly increased in the severe dengue compared to acute dengue fever."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:12:48.537137Z",
     "start_time": "2020-02-24T19:12:48.117216Z"
    }
   },
   "outputs": [],
   "source": [
    "# Pull out genes to label\n",
    "labels = [\"IFITM3\",\"IFI6\",\"IFIT3\",\"IFITM1\", \"ISG15\", \"ISG20\", \"STAT2\", \"IRF7\", \"OAS3\" ]\n",
    "labels_df = cd4_df.loc[cd4_df[\"names\"].isin( labels )]\n",
    "labels_df[\"logfoldchanges\"] = labels_df[\"logfoldchanges\"] + 0.1\n",
    "labels_df[\"neglogpval_adj\"] = labels_df[\"neglogpval_adj\"] + 0.5\n",
    "\n",
    "# Perform the plotting\n",
    "plt.figure( dpi=100, figsize=(5,5) )\n",
    "ax=plt.subplot()\n",
    "sns.scatterplot( cd4_df[\"logfoldchanges\"], cd4_df[\"neglogpval_adj\"], ax=ax )\n",
    "plt.xlabel( \"logFC\")\n",
    "plt.ylabel( \"-logP-value\")\n",
    "plt.title( \"CD4-C2-Activated: Severe vs. Acute Dengue\")\n",
    "for i in labels_df.iterrows():\n",
    "    ax.annotate( i[1][\"names\"], (i[1][\"logfoldchanges\"], i[1][\"neglogpval_adj\"]) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We generally see the same genes that were upregulated in the CD8 cells upregulated in the activated CD4 cells. This might provide evidence for the cytokine storm but, given that the same genes were upregulated in both clusters, it might just be a result of a batch effect between the samples. It could also be argued that the inteferon response genes are upregulated because of the increased viremia associated with severe dengue. More evidence is obviously needed in either case. Explore the `cd8_df` and `cd4_df` dataframes to look at what other genes are upregulated.\n",
    "\n",
    "We can additionally look at the expression of some of these genes in the three cohorts in order to determine whether they would be good biomarkers for severe dengue fever. Despite the high number of dengue infections across the world, there are no biomarkers to predict disease progression. Looking at expression in this manner is the first step in identifying which immune factors might be a good targets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-02-24T19:12:49.204886Z",
     "start_time": "2020-02-24T19:12:48.538727Z"
    }
   },
   "outputs": [],
   "source": [
    "sc.pl.violin( scdata, keys=[\"IFITM3\", \"IFIT3\", \"IFI6\"], use_raw=True, groupby=\"Disease_Status\", rotation=90 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Among the top hits, expression of *IFIT3* and *IFI6* in T cells appear to be the most suitable for use as biomarkers. Both differentiate severe dengue fever from actue cases and healthy individuals.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Summary:\n",
    "In this tutorial we went throught the typical steps of a scRNA-seq analysis workflow. We started with general preprocessing of the data which included cell and gene quality control, normalization, batch correction, selection of highly variable genes, and visualization. In these steps we explored the data and filtered it so as to conservatively remove low quality data and optimize downstream analysis results. In downstream analyses we selected T cells through clustering, and characterized T cell subsets through subclustering and differentialy gene analysis. We then used these characterized T cell subsets to inform our analysis of dengue virus infection using compositional analysis and differential gene expression. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}