Tutorial_Template_Aging_Assignment.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tutorial - Introduction to scRNA-seq analysis using Scanpy\n",
    "This Tutorial was guided heavily by Malte Luecken's \"Best Practices\" Notebook. (PostDoc from TheisLab) <br>\n",
    "The original notebook can be found at https://github.com/theislab/single-cell-tutorial/blob/master/latest_notebook/Case-study_Mouse-intestinal-epithelium_1906.ipynb\n",
    "\n",
    "It is highly recommended to look at the Scanpy Documentation to see a more detailed explanation on the function and other parameters one can specifiy: <br>\n",
    "https://scanpy.readthedocs.io/en/stable/  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Webtool to Longitudinal single cell transcriptomics:\n",
    "https://theislab.github.io/LungInjuryRegeneration/\n",
    "### Webtool to IPF Cell Atlas:\n",
    "http://www.ipfcellatlas.com/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Import Libraries used in this script\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scanpy as sc\n",
    "import pandas as pd\n",
    "import seaborn as sb\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "## Define a nice colour map for gene expression\n",
    "from matplotlib import colors\n",
    "gray_red = colors.LinearSegmentedColormap.from_list(\"\", [\"lightgray\", \"red\", \"darkred\"], N = 128)\n",
    "\n",
    "sc.settings.figdir = \"Plots/\"\n",
    "sc.set_figure_params(vector_friendly = True)\n",
    "plt.rcParams.update({\"font.size\": 14})\n",
    "sb.set_style(\"ticks\")\n",
    "\n",
    "sc.logging.print_version_and_date()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read in Data File\n",
    "The raw count matrix is stored as a .h5ad file, which can be read in directly using Scanpy’s read() function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata = ## read it in ##\n",
    "adata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Show the first 10 gene names (stored in adata.var_names) ###\n",
    "print(<TO_INSERT>)\n",
    "\n",
    "### Show the first 10 cell ids (stored in adata.obs_names) ###\n",
    "print(<TO_INSERT>)\n",
    "\n",
    "## Show the expression values of the first 10 genes in the first 5 cells\n",
    "## adata objects store the expression matrix in .X layer (cells (rows) x genes (columns))\n",
    "adata.X[0:5, 0:10].todense()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Show how many cells you have per age group, which is stored in adata.obs.grouping\n",
    "## Use the function .value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preprocessing and Quality Control\n",
    "Data quality control can be split into cell QC and gene QC. Typical quality measures for assessing the quality of a cell include the number of molecule counts (UMIs), the number of expressed genes, and the fraction of counts that are mitochondrial. A high fraction of mitochondrial reads being picked up can indicate cell stress, as there is a low proportion of nuclear mRNA in the cell. It should be noted that high mitochondrial RNA fractions can also be biological signals indicating elevated respiration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Quality control - plot QC Metrics\n",
    "plt.rcParams['figure.figsize'] = (8, 6)\n",
    "\n",
    "## use scanpy’s violin plot to visualize the quality of the cells per \"identifier\" ##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Data quality summary plots\n",
    "sc.pl.scatter(adata, \"n_counts\", \"n_genes\", color = \"identifier\")\n",
    "\n",
    "## Zoom in to range\n",
    "sc.pl.scatter(adata[adata.obs[\"n_counts\"] < 8000], \"n_counts\", \"n_genes\", color = \"percent_mito\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can assess whether there are cells with unexpected summary statistics. For example, there is a cloud of points with many counts, but few genes. This cloud of cells may indicate empty droplets containing ambient RNA, or simply dying cells.\n",
    "\n",
    "Cells with lower counts and genes tend to have a higher fraction of mitochondrial counts. These cells are likely under stress or are dying. When apoptotic cells are sequenced, there is less mRNA to be captured in the nucleus, and therefore fewer counts overall, and thus a higher fraction of counts fall upon mitochondrial RNA. If cells with high mitochondrial activity were found at higher counts/genes per cell, this would indicate biologically relevant mitochondrial activity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Thresholding decision: counts\n",
    "sb.distplot(adata.obs['n_counts'], kde = False)\n",
    "plt.show()\n",
    "\n",
    "## Zoom into lower end of the distribution\n",
    "sb.distplot(adata.obs['n_counts'][adata.obs['n_counts'] < 5000], kde = False, bins=60)\n",
    "plt.axvline(<TO_INSERT>, color = \"blue\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Thresholding decision: genes\n",
    "sb.distplot(adata.obs['n_genes'], kde=False, bins=60)\n",
    "plt.show()\n",
    "\n",
    "## Zoom into lower end of the distribution\n",
    "sb.distplot(adata.obs['n_genes'][adata.obs['n_genes'] < 3000], kde = False, bins = 60)\n",
    "plt.axvline(<TO_INSERT>, color = \"blue\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Actually do the Filtering\n",
    "print('Total number of cells: {:d}'.format(adata.n_obs))\n",
    "\n",
    "sc.pp.filter_cells(adata, max_counts = <TO_INSERT>)\n",
    "print('Number of cells after max count filter: {:d}'.format(adata.n_obs))\n",
    "\n",
    "sc.pp.filter_cells(adata, min_counts = <TO_INSERT>)\n",
    "sc.pp.filter_cells(adata, min_genes = <TO_INSERT>)\n",
    "print('Number of cells after min count filter: {:d}'.format(adata.n_obs))\n",
    "\n",
    "sc.pp.filter_genes(adata, min_cells = 3)\n",
    "print('Number of cells after gene filter: {:d}'.format(adata.n_obs))\n",
    "\n",
    "adata = adata[adata.obs[\"percent.mito\"] < 0.1]\n",
    "print('Number of cells after MT filter: {:d}'.format(adata.n_obs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The filtering is performed based on the thresholds we identified from the QC plots. Genes are also filtered if they are not detected in at least 3 cells. This reduces the dimensions of the matrix by removing 0 count genes and genes which are not sufficiently informative of the dataset. <br>\n",
    "In general it is a good idea to be permissive in the early filtering steps, and then come back to filter out more stringently when a clear picture is available of what would be filtered out. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quality Plots after Filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Re-evaluate your quality using similar plotting methods as above ### "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Summary Statistics after Filtering\n",
    "info = pd.DataFrame(data = adata.obs[\"identifier\"].cat.categories, columns = [\"identifier\"])\n",
    "info[\"n_counts\"] = adata.obs.groupby([\"identifier\"])[\"n_counts\"].median().values\n",
    "info[\"n_genes\"] = adata.obs.groupby([\"identifier\"])[\"n_genes\"].median().values\n",
    "info[\"percent_mito\"] = adata.obs.groupby([\"identifier\"])[\"percent_mito\"].mean().values\n",
    "info[\"n_cells\"] = adata.obs.groupby([\"identifier\"])[\"n_genes\"].size().values\n",
    "info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normalization and log Transformation\n",
    "Up to this point the data is only available as a count matrix. Counts are representative of molecules that were captured in the scRNA-seq experiment. As not all mRNA molecules in a cell are captured, there is a variability in the total number of counts detected between cells that results from both the number of molecules that were in the cells, and the sampling. As we cannot assume that all cells contain an equal number of molecules (cell sizes can differ substantially), we have to estimate the number of molecules that were initially in the cells. We keep this data in the '.X' part of the AnnData object as it will be used to visualize gene expression and perform statistical tests such as computing marker genes for clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Keep the count data in a counts layer \n",
    "adata.layers[\"counts\"] = adata.X.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.layers[\"counts\"][0:5, 0:10].todense()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Normalization (divide by total UMI counts per cell)\n",
    "### look at scanpy’s preprocessing methods for normalization ### "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Log Transformation\n",
    "### look at scanpy’s preprocessing methods for Log Transformation ### "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## show number of cells x number of genes\n",
    "adata.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract highly variable Genes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We extract highly variable genes (HVGs) to further reduce the dimensionality of the dataset and include only the most informative genes. Genes that vary substantially across the dataset are informative of the underlying biological variation in the data. As we only want to capture biological variation in these genes, we select highly variable genes after normalization and batch correction. HVGs are used for clustering, trajectory inference, and dimensionality reduction/visualization, while the full data set is used for computing marker genes, differential testing and visualizing expression values on the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.highly_variable_genes(adata, flavor = \"seurat\", n_top_genes = 4000, batch_key = \"identifier\")\n",
    "print('\\n','Number of highly variable genes: {:d}'.format(np.sum(adata.var['highly_variable'])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Find a visualization to show characeristics of the chosen genes ###"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Highly variable gene information is stored automatically in the adata.var['highly_variable'] field. The dataset now contains:\n",
    "\n",
    "- log-normalized data in adata.X\n",
    "- highly variable gene annotations in adata.var['highly_variable']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dimensionality Reduction using Principal Component Analysis\n",
    "Visualizing scRNA-seq data is the process of projecting a high-dimensional matrix of cells and genes into a few coordinates such that every cell is meaningfully represented in a two-dimensional graph. Tt is a good idea to look at several visualizations and decide which visualization best represents the aspect of the data that is being investigated. <br>\n",
    "Overall t-SNE visualizations have been very popular in the community, however the recent UMAP algorithm has been shown to better represent the topology of the data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Perform a Principal Component Analysis using scanpy’s preprocessing methods ###\n",
    "### Remember to specify to use only the previously selected highly variable genes for this step ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.pca_variance_ratio(adata, n_pcs = 40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Calculate nearest neighbourhood graph and find fitting values for the parameter n_neighbors, n_pcs ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Generate a 2D embedding using UMAP and plot the result ###\n",
    "plt.rcParams['figure.figsize'] = (6, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## You can see the annotation of each cell (can be used to colour the cells) as the columns of\n",
    "adata.obs.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save the Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Find out how to save the plot as PDF. Hint: look at the save parameter of sc.pl.umap()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clustering\n",
    "Clustering is a central component of the scRNA-seq analysis pipeline. To understand the data, we must identify cell types and states present. The first step of doing so is clustering. This modularity optimization is performed on the k-nearest-neighbour graph of cells in scRNA-seq analysis.  <br>\n",
    "By changing the resolution parameter one can influence the number of clusters the algorithm finds. Investigating several resolutions allows us to select a clustering that appears to capture the main clusters in the visualization and can provide a good baseline for further subclustering of the data to identify more specific substructures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Generate an unsupervised Clustering using scanpy’s tl.leiden() function ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## You can increase the plot by setting following Paramter in the funtion\n",
    "ax = plt.figure(figsize = (5, 4), dpi = 120).gca()\n",
    "\n",
    "## One can directly plot the labels on top of the cells by setting\n",
    "legend_loc = \"on data\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can also visualize expression of certain known marker genes by specifying the colour Parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color = [\"Ear2\", \"Sftpd\", \"Foxj1\"], cmap = gray_red, size = 20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find Marker Genes and Annotate Clusters\n",
    "To annotate the clusters we obtained, we find genes that are up-regulated in the cluster compared to all other clusters (marker genes). The test is automatically performed on the .raw data set, which is uncorrected and contains all genes. All genes are taken into account, as any gene may be an informative marker."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Generate a markers table using scanpy’s rank_genes_groups() and specify the group info ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Plot the first 5 genes per cluster\n",
    "sc.pl.rank_genes_groups_dotplot(adata, n_genes = 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To further identify the clusters in ones data set, one can look at the overlap with a list of known marker genes.\n",
    "In practice marker gene sets can be obtained from public databases such as Linnarson's mouse brain atlas, various Human Cell Atlas datasets, and other published reference atlases. It should be noted that marker genes may not always overlap as expected given that atlases tend to be investigated under wild-type conditions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pick one or two cluster for which you could identify the cell type and annotate them in object\n",
    "As an example, here we are showing the expression of 3 well known marker genes for cell types of the lung.  \n",
    "Feel free to use other marker genes that you have in mind, or to check out the marker genes saved in Plot/dotplot_cell_type_marker.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.violin(adata, keys = [\"Sftpd\"], groupby = \"leiden_1\", size = 1.5)    ## AT2 cell marker\n",
    "sc.pl.violin(adata, keys = [\"Foxj1\"], groupby = \"leiden_1\", size = 1.5)    ## Ciliated cell marker\n",
    "sc.pl.violin(adata, keys = [\"Ear2\"], groupby = \"leiden_1\", size = 1.5)     ## Macrophage marker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.dotplot(adata, var_names = [\"Ear2\", \"Sftpd\", \"Foxj1\"], groupby = \"leiden_1\", dendrogram = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## With this code snippet, you can highlight certain cluster of interest in the UMAP\n",
    "sc.pl.umap(adata, color = \"leiden_1\", legend_loc = \"on data\", \n",
    "           groups = [\"1\", \"17\", \"8\", \"14\", \"20\", \"6\", \"0\", \"10\", \"4\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## You can pick certain cluster and assign them a cell type as following\n",
    "celltype_map = {\"1\": \"Macrophages\", \"4\": \"Ciliated cells\", \"17\": \"AT2 cells\", \"8\": \"AT2 cells\",\n",
    "                \"14\": \"AT2 cells\", \"20\": \"AT2 cells\", \"6\": \"AT2 cells\", \"0\": \"AT2 cells\", \"10\": \"AT2 cells\",\n",
    "                \"17\": \"Doublets\"}\n",
    "\n",
    "adata.obs[\"cell_type\"] = [celltype_map[leiden] if leiden in celltype_map.keys() else leiden\n",
    "                          for leiden in adata.obs.leiden]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color = \"cell_type\", legend_loc = \"on data\", ax = plt.figure(figsize=(6, 5), dpi = 100).gca())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pick one Example Cluster for Differential Gene Expression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Create a new object, using only cells assigned to\n",
    "cell_type = <TO_INSERT>\n",
    "sub = adata[adata.obs[\"cell_type\"] == cell_type].copy()\n",
    "sc.pl.umap(sub, color = [\"cell_type\", \"grouping\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After subsetting the data, it makes sense to re-caculate the PCs and knn graph, as the neighbourhood of the individual cells has likely changed. This also helps enhance certain features which were over-shadowed in the full embedding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### recalculate PCA, k nearest neighbor graph and the UMAP ###\n",
    "### Keep in mind to use the right object in the functions :) ###"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perform Differential Gene Expression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Use tl.rank_genes_groups to find genes upregulated in either young or old mice ###\n",
    "## Examine the results e.g. by sc.pl.rank_genes_groups_dotplot ###"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Potential Bonus: PAGA\n",
    "By quantifying the connectivity of partitions (groups, clusters) of the single-cell graph, partition-based graph abstraction (PAGA) generates a much simpler abstracted graph of partitions, in which edge weights represent confidence in the presence of connections. By tresholding this confidence in paga(), a much simpler representation of the manifold data is obtained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Perform PAGA on the whole object ###"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### You can save the object and all annotations calculated so far with the following command"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.write(\"Data/Aging_processed.h5ad\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}