module-5-scrnaseq/scRNAseq_Winter_2022.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Single Cell RNA-seq Analysis\n",
    "## scRNA-seq experiments of day 21 reprogramming intermediates cultured under fibroblast condition, naive pluripotent and trophoblast stem cell conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#something like:\n",
    "from IPython.display import Image\n",
    "Image(\"Title.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This dataset was downloaded from: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE150637"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#something like:\n",
    "from IPython.display import Image\n",
    "Image(\"Experimental_Design.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Import Python packages/libraries to use in analysis¶"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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",
    "from matplotlib import rcParams\n",
    "from matplotlib import colors\n",
    "import seaborn as sns\n",
    "import warnings;\n",
    "warnings.filterwarnings('ignore');\n",
    "from gprofiler import GProfiler\n",
    "sc.settings.verbosity = 3 # verbosity: errors (0), warnings (1), info (2), hints (3)\n",
    "sc.settings.set_figure_params(dpi=170, color_map='viridis')  # low dpi (dots per inch) yields small inline figures\n",
    "sc.logging.print_header()\n",
    "results_file = './liu2020.h5ad'\n",
    "results_file_denoised = './liu2020_deno.h5ad'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Note that we need to: pip install fa2\n",
    "- Open a terminal window on datahub\n",
    "- Activate the conda environment for this module: `conda activate scanpy_2021`\n",
    "- Install fa2 using: `pip install fa2`\n",
    "- Deactivate the environment: `conda deactivate`\n",
    "- Return to the notebook!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we import a package as something else, we are creating a shortcut to easily call the package. For example, import pandas as pd will allow us to type pd when calling a pandas function rather than pandas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Starting from a 10X Dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we opened jupyter lab/notebook from the directory that we made, we only need to enter the path and filename starting from this point. We are going to use Scanpy's \"read_10x_mtx\" function because we are reading in a 10X mtx file to analyze. We could also start with a matrix generated from FeatureCounts or DESeq2 or any other program using the python package pandas or use a .h5ad file from cell ranger (10X).\n",
    "​\n",
    "\n",
    "Here, we are reading in sparse count matrices. If we were going to perform batch correction we would need to convert these to dense representation with the .toarray() function. We create dense matrices as our batch correction method outputs a dense expression matrix, and the data transfer between R and python is currently limited to dense matrices. When sparse batch correction methods are available, and rpy2 is extended to sparse matrices, it is more memory-efficient to keep the data in a sparse format.\n",
    "​\n",
    "\n",
    "It should also be noted that the conventions for storing single cell data differ between R (Seurat, or Scater) and python platforms (Scanpy). Scanpy expects the data to be stored in the format (cells x genes), while R platforms expect the transpose of that. As data is typically stored in the format (genes x cells) as in GEO database, we must transpose the data before using it. However, because we are using scanpy's built in 10x data loading, it will transpose the data for us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_TSC = sc.read_10x_mtx('./TSC/hg19/',  # the directory with the `.mtx` file\n",
    "    var_names='gene_symbols',                  # use gene symbols for the variable names (variables-axis index)\n",
    "    cache=True)                                # write a cache file for faster subsequent reading"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The command \"read_10x_mtx\" assumes that the data in question is named matrix.mtx and that there are \"features\" and \"barcodes\" .tsv files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_TSC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above shows that we have about 4700 cells and almost 33,000 genes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can't have duplicate gene names so we need to make sure that they are unique"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_TSC.var_names_make_unique()  # this is unnecessary if using 'gene_ids'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_TSC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_HDF1 = sc.read_10x_mtx('./HDF1/hg19/',  # the directory with the `.mtx` file\n",
    "    var_names='gene_symbols',                  # use gene symbols for the variable names (variables-axis index)\n",
    "    cache=True)                                # write a cache file for faster subsequent reading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_HDF1.var_names_make_unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_HDF1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_naive = sc.read_10x_mtx('./Naive/hg19/',  # the directory with the `.mtx` file\n",
    "    var_names='gene_symbols',                  # use gene symbols for the variable names (variables-axis index)\n",
    "    cache=True)                                # write a cache file for faster subsequent reading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_naive.var_names_make_unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_naive"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now I am going to label all of the observations from each of the separate data sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_TSC.obs['cell_type'] = 'TSC'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_HDF1.obs['cell_type'] = 'HDF1'\n",
    "adata_naive.obs['cell_type'] = 'Naive'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_TSC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now I am going to concetenate the datasets together and start the preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata = adata_TSC.concatenate([adata_HDF1, adata_naive])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 Pre-processing and Visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 Quality Control"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.\n",
    "\n",
    "Alternatively, preprocessing can be done somewhere else (such as in R) and then loaded into scanpy. If this is your situation, skip this portion. You can also freeze the raw data at any point if you want to come back to it. Set the .raw attribute of AnnData object to the logarithmized raw gene expression for later use in differential testing and visualizations of gene expression. This simply freezes the state of the AnnData object returned by sc.pp.log1p.\n",
    "\n",
    "Log transforming the data in one simple way to \"variance-stabilize\" your data to adjust the data so that the influence of gene expression magnitude on gene variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Quality control - calculate QC covariates\n",
    "adata.obs['n_counts'] = sum(adata.X.sum(1).tolist(),[])\n",
    "adata.obs['log_counts'] = np.log(adata.obs['n_counts'])\n",
    "adata.obs['n_genes'] = sum((adata.X > 0).sum(1).tolist(),[])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Citing from “Simple Single Cell” workflows (Lun, McCarthy & Marioni, 2017):\n",
    "\n",
    "High proportions of mitochondrial reads are indicative of poor-quality cells (Islam et al. 2014; Ilicic et al. 2016), possibly because of loss of cytoplasmic RNA from perforated cells. The reasoning is that mitochondria are larger than individual transcript molecules and less likely to escape through tears in the cell membrane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mito_genes = adata.var_names.str.startswith('MT-')\n",
    "# for each cell compute fraction of counts in mito genes vs. all genes\n",
    "# the `.A1` is only necessary as X is sparse (to transform to a dense array after summing)\n",
    "adata.obs['percent_mito'] = np.sum(adata[:, mito_genes].X, axis=1).A1 / np.sum(adata.X, axis=1).A1\n",
    "# add the total counts per cell as observations-annotation to adata\n",
    "adata.obs['n_counts'] = adata.X.sum(axis=1).A1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obs.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scanpy has a \"cheat\" command called \"pp.calculate_qc_metrics()\" that we could also use instead of doing all the calculations we did above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#first we need to make a copy so we don't write anything in our Anndata\n",
    "adata_qcTest = adata.copy()\n",
    "sc.pp.calculate_qc_metrics(adata_qcTest, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#check to see we got the same numbers as our \"manual\" calculations above:\n",
    "adata_qcTest.obs.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets use a violin plot to take a look at the quality of our data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.violin(adata, ['n_genes', 'n_counts'],jitter=0.4)\n",
    "sc.pl.violin(adata, ['percent_mito'],jitter=0.4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plots show the counts per cell and the fraction of mitochondrial reads (MT frac) for most cells are still far below 20-25%, which is the typical filtering thresholds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Data quality summary plots\n",
    "p1 = sc.pl.scatter(adata, 'n_counts', 'n_genes', color='percent_mito')\n",
    "p2 = sc.pl.scatter(adata[adata.obs['n_counts']<10000], 'n_counts', 'n_genes', color='percent_mito')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By looking at plots of the number of genes versus the number of counts with MT fraction information, we can assess whether there are cells with unexpected summary statistics. It is important here to look at these statistics jointly. For example, there is a cloud of points with many counts, but few genes. Our first instinct would be to filter out these as \"dying\" outliers, however this should only be done if they show high MT fraction. In this case we should probably filter out the cells with very few genes as these may be difficult to annotate later. This will be true for the initial cellular density between 1000-4000 counts and < ~500 genes.\n",
    "\n",
    "It is typical for the main cloud of data points, that 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",
    "p3 = sns.distplot(adata.obs['n_counts'], kde=False)\n",
    "plt.show()\n",
    "\n",
    "p4 = sns.distplot(adata.obs['n_counts'][adata.obs['n_counts']<10000], kde=False, bins=60)\n",
    "plt.show()\n",
    "\n",
    "p5 = sns.distplot(adata.obs['n_counts'][adata.obs['n_counts']>25000], kde=False, bins=60)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Histograms of the number of counts per cell show that there is a groups of cells with fewer than 1,500 counts. These are likely uninformative given the overall distribution of counts. This may be cellular debris found in droplets.\n",
    "\n",
    "On the upper end of counts, we see just a few cells with counts over 50,000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Thresholding decision: genes\n",
    "p6 = sns.distplot(adata.obs['n_genes'], kde=False, bins=60)\n",
    "plt.show()\n",
    "\n",
    "p7 = sns.distplot(adata.obs['n_genes'][adata.obs['n_genes']<2500], kde=False, bins=60)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like our data has a small population of cells with low gene counts. Given these plots, and the plot of genes vs counts above, we decide to filter out cells with fewer than 1500 genes expressed. Below this we are likely to find dying cells or empty droplets with ambient RNA. Looking above at the joint plots, we see that we filter out the main density of low gene cells with this threshold\n",
    "\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. This is available after visualization/clustering. For demonstration purposes we stick to a simple (and slightly more stringent) filtering here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter cells according to identified QC thresholds:\n",
    "print('Total number of cells: {:d}'.format(adata.n_obs))\n",
    "\n",
    "sc.pp.filter_cells(adata, min_counts = 1500)\n",
    "print('Number of cells after min count filter: {:d}'.format(adata.n_obs))\n",
    "\n",
    "sc.pp.filter_cells(adata, max_counts = 50000)\n",
    "print('Number of cells after max count filter: {:d}'.format(adata.n_obs))\n",
    "\n",
    "adata = adata[adata.obs['percent_mito'] < 0.2]\n",
    "print('Number of cells after MT filter: {:d}'.format(adata.n_obs))\n",
    "\n",
    "sc.pp.filter_cells(adata, min_genes = 1500)\n",
    "print('Number of cells after gene filter: {:d}'.format(adata.n_obs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Filter genes:\n",
    "print('Total number of genes: {:d}'.format(adata.n_vars))\n",
    "\n",
    "# Min 15 cells - filters out 0 count genes\n",
    "sc.pp.filter_genes(adata, min_cells=30)\n",
    "print('Number of genes after cell filter: {:d}'.format(adata.n_vars))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Genes are also filtered if they are not detected in at least 30 cells. This reduces the dimensions of the matrix by removing 0 count genes and genes which are not sufficiently informative of the dataset. 30 is not a magic number and should be adjusted down or up depending on how many cells are in your dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 Normalization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each count in a count matrix represents the successful capture, reverse transcription and sequencing of a molecule of cellular mRNA. Count depths for identical cells can differ due to the variability inherent in each of these steps. Thus, when gene expression is compared between cells based on count data, any difference may have arisen solely due to sampling effects. Normalization addresses this issue by e.g. scaling count data to obtain correct relative gene expression abundances between cells.\n",
    "\n",
    "Scanpy has a builtin command (scanpy.pp.normalize_total) which normalizes the counts per cell. If we choose target_sum=1e6, this is CPM normalization. CPM normalization is the most commonly used normalization protocol and is also referred to as “counts per million”. This protocol comes from bulk expression analysis and normalizes count data using a so-called size factor proportional to the count depth per cell. Variations of this method scale the size factors with different factors of 10, or by the median count depth per cell in the dataset. CPM normalization assumes that all cells in the dataset initially contained an equal number of mRNA molecules and count depth differences arise only due to sampling. If target_sum=None, after normalization, each observation (cell) has a total count equal to the median of total counts for observations (cells) before normalization \n",
    "\n",
    "If exclude_highly_expressed=True, very highly expressed genes are excluded from the computation of the normalization factor (size factor) for each cell. This is meaningful as these can strongly influence the resulting normalized values for all other genes [Weinreb17]\n",
    "\n",
    "We will use this normalization for our data today because we have limited cellular variability in our data and due to the time and effort it takes to set up SCRAN, but it is probably not the best method to use."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Important: Before normalizing the data, we ensure that a copy of the raw count data is kept in a separate AnnData object. This allows us to use methods downstream that require this data as input.**"
   ]
  },
  {
   "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": [
    "#Normalize the data and log transform it\n",
    "sc.pp.normalize_per_cell(adata)\n",
    "sc.pp.log1p(adata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Store the full data set in 'raw' as log-normalised data for statistical testing\n",
    "adata.raw = adata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The count data has been normalized and log-transformed with an offset of 1. The latter is performed to normalize the data distributions. The offset of 1 ensures that zero counts map to zeros. We keep this data in the '.raw' 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": "markdown",
   "metadata": {},
   "source": [
    "## 2.3 Batch Correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Batch correction is performed to adjust for batch effects from the different samples that you might load together. Our dataset contains three different conditions applied to the same starting cells that are then being combined together so we want to correct for any batch effects. As the batch effect from samples are often overlapping, correcting for this batch effect will also partially regress out differences between regions. We allow for this to optimally cluster the data. This approach can also be helpful to find differentiation trajectories, but we revert back to non-batch-corrected data for differential testing and computing marker genes.\n",
    "\n",
    "Regression batch correction essentially assumes that each batch must have the same average gene expression across all cells and removes the difference (this can over-correct your data).\n",
    "\n",
    "Note that ComBat batch correction requires a dense matrix format as input (which is already the case in this example).\n",
    "\n",
    "If we wanted to batch correct we would perform the below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ComBat batch correction\n",
    "sc.pp.combat(adata, key='batch')\n",
    "#where \"batch\" is what want to correct for. We set this up when we load the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From Li et al. 2020 \"Although the Regress_Out and ComBat algorithms are both based on linear regression, we noticed the superiority of the ComBat method compared to the Regress_Out method in terms of the batch-correction performance and the computing efficiency (Figs. 1, 2, S9). The Regress_Out algorithm uses a general linear model (GLM) to regress out unwanted sources of variation in the expression matrix. The coefficient for each batch block is estimated by fitting the GLM and is then set to zero to remove the corresponding batch effect. The ComBat algorithm uses the same strategy but performs an empirical Bayes (EB) framework to adjust the expression matrix for batch correction. The EB framework is usually designed to shrink gene variances and remove the inferred batch effect by fitting a Bayesian model. ComBat first standardizes all gene expression values across cells and fits the data into a standard distribution Bayesian model. Then, it can infer the batch effects using the estimated model and adjust the gene variances (Johnson et al. 2007a; Leek 2014). Due to the empirical Bayesian shrinkage of the blocking coefficient estimates, the ComBat algorithm is quite robust when processing diverse batches, which is also reflected in our analysis in both the MCA and TM data sets (Figs. 1, 2).\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note ComBat batch correction can produce negative expression values. One can either set all negative values to 0 or force zeros pre-batch-correction to remain zero post-batch-correction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes batch correction isn't enough...\n",
    "There are also several data integration tools that you can and should use to integrate single cell data from multiple experiments such as bbknn (https://doi.org/10.1093/bioinformatics/btz625) or Harmony (Nature Methods <https://doi.org/10.1038/s41592-019-0619-0>). Scanpy has the \"ingest\" (scanpy.tl.ingest) built into the API but also has bbknn, Harmony, and mutual nearest neighbor in the external API (pp.bbknn(adata[, batch_key, approx, metric, …]) and (pp.harmony_integrate(adata, key[, basis, …]) and (pp.mnn_correct(*datas[, var_index, …])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.4 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, cell cycle scoring, and visualizing expression values on the data.**\n",
    "\n",
    "Keep in mind that any given tissue or cell type expresses 11,000 to 13,000 genes, of which 3,000 to 5,000 have a cell type expression pattern, whereas the remaining genes are ubiquitously expressed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTE: The below command expects log data (except when the flavor='seurat_v3') is used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=5, min_disp=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.highly_variable_genes(adata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plots show how the data was normalized to select highly variable genes irrespective of the mean expression of the genes unless  you set a maximum. This is achieved by using the index of dispersion which divides by mean expression, and subsequently binning the data by mean expression and selecting the most variable genes within each bin."
   ]
  },
  {
   "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",
    "- a 'counts' layer with count data\n",
    "- log-normalized data in adata.raw\n",
    "- highly variable gene annotations in adata.var['highly_variable']\n",
    "\n",
    "The HVG labels will be used to subselect genes for clustering and trajectory analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we could use the standard technique for the extraction of highly variable genes from the 10X genomics preprocessing software CellRanger. Typically between 1000 and 5000 genes are selected. Here, we extract the top 2000 most variable genes for further processing. If particular genes of importance are known, one could assess how many highly variable genes are necessary to include all, or the majority, of these"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.highly_variable_genes(adata, flavor='cell_ranger', n_top_genes=2000)\n",
    "print('\\n','Number of highly variable genes: {:d}'.format(np.sum(adata.var['highly_variable'])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To take a look at how many genes we ended up with using the filter above we can do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.var.highly_variable.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And to see how many genes are not HVG:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(~adata.var.highly_variable).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.5 Dimenstionality Reduction Methods: Visualization and Summarization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are two main objectives of dimensionality reduction methods: visualization and summarization. Visualization is the attempt to optimally describe the dataset in two or three dimensions. These reduced dimensions are used as coordinates on a scatter plot to obtain a visual representation of the data. Summarization does not prescribe the number of output components. Instead, higher components become less important for describing the variability present in the data. Summarization techniques can be used to reduce the data to its essential components by finding the inherent dimensionality of the data, and are thus helpful for downstream analysis.\n",
    "\n",
    "\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. However, the visualization of scRNA-seq data is an active area of research and each method defines 'meaningful' in its own way. Thus, it is a good idea to look at several visualizations and decide which visualization best represents the aspect of the data that is being investigated.\n",
    "\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. One criticism is that they neglect the local density of data points in the original space, often resulting in misleading visualiza-\n",
    "tions where densely populated subsets of cells are given more visual space than warranted by their transcriptional diversity\n",
    "in the dataset.\n",
    "\n",
    "Note that we do not scale the genes to have zero mean and unit variance. A lack of rescaling is equivalent to giving genes with a higher mean expression a higher weight in dimensionality reduction (despite correcting for mean offsets in PCA, due to the mean-variance relationship). We argue that this weighting based on mean expression being a biologically relevant signal. However, rescaling HVG expression is also common, and the number of publications that use this approach suggests that scaling is at least not detrimental to downstream scRNA-seq analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# However, there is some debate..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "Image(\"liorPachter_2021.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what Lior says on twitter:\n",
    "\"It's time to stop making t-SNE & UMAP plots. In a new preprint w/ Tara Chari we show that while they display some correlation with the underlying high-dimension data, they don't preserve local or global structure & are misleading. They're also arbitrary.Threadhttps://biorxiv.org/content/10.1101/2021.08.25.457696v1\"\n",
    "\n",
    "\n",
    "and...\n",
    "\n",
    "\n",
    "\"Ok.. but.. maybe t-SNE & UMAP (or your favorite 2D viz) aren't perfect, but they are \"canonical\" and not arbitrary. Nope. They're just art.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the visualizations\n",
    "sc.pp.pca(adata, n_comps=50, use_highly_variable=True, svd_solver='arpack')\n",
    "sc.pp.neighbors(adata)\n",
    "\n",
    "sc.tl.tsne(adata, n_jobs=12) #Note n_jobs works for MulticoreTSNE, but not regular implementation)\n",
    "sc.tl.umap(adata)\n",
    "sc.tl.diffmap(adata)\n",
    "sc.tl.draw_graph(adata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.pca_scatter(adata, color='n_counts')\n",
    "sc.pl.tsne(adata, color='n_counts')\n",
    "sc.pl.umap(adata, color='n_counts')\n",
    "sc.pl.diffmap(adata, color='n_counts', components=['1,2','1,3'])\n",
    "sc.pl.draw_graph(adata, color='n_counts')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PCA:\n",
    "- Unsurprisingly, the first principle component captures variation in count depth between cells, and is thus only marginally informative\n",
    "- The plot usually will not show the expected clustering of the data in two dimensions unless your data is not very heterogenous as ours is\n",
    "\n",
    "t-SNE:\n",
    "- Shows several distinct clusters with clear subcluster structure\n",
    "- Connections between clusters are difficult to interpret visually\n",
    "\n",
    "UMAP:\n",
    "- Data points are spread out on the plot showing several clusters\n",
    "- Connections between clusters can be readily identified\n",
    "\n",
    "Diffusion Maps:\n",
    "- Shows connections between regions of higher density\n",
    "- Very clear trajectories are suggested, but clusters are less clear\n",
    "- Each diffusion component extracts heterogeneity in a different part of the data\n",
    "\n",
    "Graph:\n",
    "- Does not show much new information when compared to the others but often will show clear connections between the clusters and hint at stem cell clusters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The strengths and weaknesses of the visualizations might not be readily identified in the above plots. However, generally t-SNE exaggerates differences (but sacrafices accurate global distance representation in favor of local distance), diffusion maps exaggerate transitions. Overall UMAP and force-directed graph drawings show the best compromise of the two aspects, however UMAP is much faster to compute. UMAP has furthermore been shown to more accurately display the structure in the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way to make the above run faster is to use the elbow method to figure out how many PC's to use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.pca_variance_ratio(adata, log=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before and after batch correction or after dimensionality reduction it is a good idea to take a look at how your different batches clustered"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['n_counts', 'batch'], alpha=0.5, size=25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.draw_graph(adata, color='batch')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.6 Cell Cycle Scoring"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far we have tried to correct for known sources of technical variation in the data (e.g. batch, count depth). There are also known sources of biological variation that can explain the data and also can be corrected for (e.g. cell cycle). This data correction can be performed by a simple linear regression against a cell cycle score as implemented in the Scanpy and Seurat platforms (Butler et al, 2018; Wolf et al, 2018) or in specialized packages with more complex mixture models such as scLVM (Buettner et al, 2015) or f-scLVM (Buettner et al, 2017). Lists of marker genes to compute cell cycle scores are obtained from the literature (Macosko et al, 2015). These methods can also be used to regress out other known biological effects such as mitochondrial gene expression, which is interpreted as an indication of cell stress. Several aspects should be considered prior to correcting data\n",
    "for biological effects. Firstly, correcting for biological covariates is not always helpful to interpret scRNA-seq data. While removing cell cycle effects can improve the inference of developmental trajectories (Buettner et al, 2015; Vento-Tormo et al, 2018), cell cycle signals can also be informative of the biology. \n",
    "\n",
    "Here, we will use a gene list from Macosko et al., Cell (2015) to score the cell cycle effect in the data and classify cells by cell cycle phase. Please note, the gene list was generated for human HeLa cells. We perform cell cycle scoring on the full batch-corrected data set in adata."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Score cell cycle and visualize the effect:\n",
    "cc_genes = pd.read_table('./Macosko_cell_cycle_genes.txt', delimiter='\\t')\n",
    "s_genes = cc_genes['S'].dropna()\n",
    "g2m_genes = cc_genes['G2.M'].dropna()\n",
    "\n",
    "s_genes_ens = adata.var_names[np.in1d(adata.var_names, s_genes)] #in1d = Test whether each element of a 1-D array is also present in a second array\n",
    "g2m_genes_ens = adata.var_names[np.in1d(adata.var_names, g2m_genes)]\n",
    "\n",
    "sc.tl.score_genes_cell_cycle(adata, s_genes=s_genes_ens, g2m_genes=g2m_genes_ens)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['S_score', 'G2M_score'], use_raw=False)\n",
    "sc.pl.umap(adata, color='phase', use_raw=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "None of our clusters show a pronounced signal and look like they are dividing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we did want to regress out a set of genes we would do something like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sc.pp.regress_out(adata, ['S_score', 'G2M_score'])\n",
    "#sc.pp.scale(adata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3 Downstream analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1 Clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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. Performing Modularity optimization by Louvain community detection on the k-nearest-neighbour graph of cells has become an established practice in scRNA-seq analysis. Thus, this is the method of choice in this tutorial as well.\n",
    "\n",
    "Here, we perform clustering at two resolutions. 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 substructure.\n",
    "\n",
    "Clustering is performed on the highly variable gene data, dimensionality reduced by PCA, and embedded into a KNN graph. (see sc.pp.pca() and sc.pp.neighbors() functions used in the visualization section.\n",
    "\n",
    "\n",
    "Modularity optimization via louvain has a stochastic element to it. This stochasticity typically does not affect the biological interpretation of the data, but can change the details of analysis scripts. Normally scanpy fixes the random seed to 0 to make scripts exactly reproducible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.louvain(adata, key_added='louvain_r1')\n",
    "sc.tl.louvain(adata, resolution=0.5, key_added='louvain_r0.5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obs['louvain_r0.5'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obs['louvain_r1'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#Visualize the clustering and how this is reflected by different technical covariates\n",
    "sc.pl.umap(adata, color=['louvain_r1', 'louvain_r0.5'], palette='Set2')\n",
    "sc.pl.umap(adata, color=['n_counts'])\n",
    "sc.pl.umap(adata, color=['batch'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At a resolution of 0.5 the broad clusters in the visualization are captured well in the data. It remanins to be seen if cluster 2 in the 0.5 resolution should be two or three clusters instead of one (as in the 1 resolution). The covariate plots show that none of the clusters have unusually high or low counts, are made up of just the batches, and we could also check for a high percent mitochondria if we wanted to determine if some of the clusters seem to represent stressed or dying cells."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To look at differences in clustering algorithms we will cluster cells using the Leiden algorithm [Traag18], an improved version of the Louvain algorithm [Blondel08]. It has been proposed for single-cell analysis by [Levine15].\n",
    "\n",
    "It might also be interesting to take a look at 0.75 resolution as well"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.louvain(adata, resolution=0.75, key_added='louvain_r0.75')\n",
    "sc.tl.leiden(adata, key_added='leiden_r1')\n",
    "sc.tl.leiden(adata, resolution=0.5, key_added='leiden_r0.5')\n",
    "sc.tl.leiden(adata, resolution=0.75, key_added='leiden_r0.75')\n",
    "sc.pl.umap(adata, color=['louvain_r1', 'leiden_r1'], palette='Set2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now I want to see what the 0.75 resolution looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['louvain_r0.75', 'leiden_r0.75'], palette='Set2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no \"right\" answer here but I think the Leiden 0.75 looks more resonable. Most likely, you will need to do several different methods and then use your biological knowledge to make the decision."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next let's look at the \"cell type\". Remember this is what we named the cells when we first uploaded them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=['cell_type', 'leiden_r0.75'], palette='Set2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.draw_graph(adata, color=['cell_type', 'leiden_r0.75'], palette='Set2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2 Marker genes & cluster annotation "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To annotate the clusters we obtained, we find genes that are up-regulated in the cluster compared to all other clusters (marker genes). This differential expression test is performed by a Welch t-test. This is the default in scanpy. 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.\n",
    "\n",
    "Also note that the results are displayed as ranked gene lists instead of genes that have a p-value below some cut-off. This is because the test treats each cell as an independent replicate leading to very low p-values. It might be useful in some cases to create pseudo-bulk counts.\n",
    "\n",
    "We can use some known marker genes to annotate the clusters and if we didn't have any information on the cells in question we could run automated annotation using scmap or garnett or the recently published MARS (https://github.com/snap-stanford/mars) as well as do a gene ontolody enrichment analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Calculate marker genes\n",
    "sc.tl.rank_genes_groups(adata, groupby='leiden_r0.75', key_added='rank_genes_r0.75')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot marker genes\n",
    "sc.pl.rank_genes_groups(adata, key='rank_genes_r0.75', groups=['0','1','2'], fontsize=12)\n",
    "sc.pl.rank_genes_groups(adata, key='rank_genes_r0.75', groups=['3','4','5'], fontsize=12)\n",
    "sc.pl.rank_genes_groups(adata, key='rank_genes_r0.75', groups=['6','7','8'], fontsize=12)\n",
    "sc.pl.rank_genes_groups(adata, key='rank_genes_r0.75', groups=['9','10','11'], fontsize=12)\n",
    "sc.pl.rank_genes_groups(adata, key='rank_genes_r0.75', groups=['12','13'], fontsize=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clusters 4, 6, and 11 all look similar. If we wanted to find the differences between these three clusters we could do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.rank_genes_groups(adata,groupby='leiden_r0.75',groups=['6'],reference='11',key_added='rank_genes_r0.75_clust6')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.rank_genes_groups(adata, key='rank_genes_r0.75_clust6', groups=['6'], n_genes=30, sharey=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make things easier to look at we can put everything in a table. Here are the top 20 genes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.DataFrame(adata.uns['rank_genes_r0.75']['names']).head(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we would like to take a look at some of the marker genes across the groups we can do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.violin(adata, ['ALPPL2', 'GATA3', 'COL1A1'], groupby='leiden_r0.75')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From here, we can use our biological knowledge to start labeling some of the clusters. However, we have a list from the Liu et al. paper (Supplementary File 3). They write in the methods: \"To identify gene signatures, we first define cluster-specific marker genes for each of the 21 snRNA clusters. For each snRNA cluster, we define marker genes as genes that have an average log 2 -transformed fold change (averaged across all 20 pairwise differential expressions) of >1.5, and we also require the genes to be differentially expressed in at least 14 of the 20 pairwise differential expressions. Hierarchical clustering was then performed on the Jaccard similarity of the marker genes (Extended Data Fig. 2f) to identify overlapping gene sets that is, the gene signatures. Overall, we identified eight gene signatures (Supplementary Table 3) that we named fibroblast (snRNA-FM1, snRNA-FM2, snRNA-FM3, snRNA-FM4); mixed (snRNA-mix); early-primed (snRNA-PR1); primed: snRNA-PR2, snRNA-PR3, snRNA-PR4); novel intermediate signature (snRNA-nic); naive (snRNA-NR1, snRNA-NR2, snRNA-NR3, snRNA-NR4); non-reprogrammed (nonReprog)1 (snRNA-RE1, snRNA-RE3, snRNA-RE4, snRNA-RE5); nonReprog2 (snRNA-RE6).\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use these gene signatures to score each of the clusters and identify them in this manner. To score a gene list, the algorithm calculates the difference of mean expression of the given list and the mean expression of reference genes. To build the reference, the function randomly chooses a bunch of genes matching the distribution of the expression of the given list."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we can calculate an overlap score between data-deriven marker genes and provided markers.\n",
    "\n",
    "Marker gene overlap scores can be quoted as overlap counts, overlap coefficients, or jaccard indices. The method returns a pandas dataframe which can be used to annotate clusters based on marker gene overlaps."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because they did something similar in Liu et al., lets give it a try and score the clusters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#First we will pull in the gene signature lists and do the same thing we did for cell cycle\n",
    "sig_genes = pd.read_table('./SupTable3.txt', delimiter='\\t')\n",
    "fibroblast = sig_genes['fibroblast'].dropna()\n",
    "mixed = sig_genes['mixed'].dropna()\n",
    "early_primed = sig_genes['early-primed'].dropna()\n",
    "primed = sig_genes['primed'].dropna()\n",
    "naive = sig_genes['naive'].dropna()\n",
    "nis = sig_genes['nis'].dropna()\n",
    "nonReprog1 = sig_genes['nonReprog1'].dropna()\n",
    "nonReprog2 = sig_genes['nonReprog2'].dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.score_genes(adata, fibroblast, score_name='fibro_score')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.draw_graph(adata, color=['cell_type', 'leiden_r0.75','fibro_score'], palette='Set2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.score_genes(adata, mixed, score_name='mix_score')\n",
    "sc.tl.score_genes(adata, early_primed, score_name='early_prime_score')\n",
    "sc.tl.score_genes(adata, primed, score_name='primed_score')\n",
    "sc.tl.score_genes(adata, naive, score_name='naive_score')\n",
    "sc.tl.score_genes(adata, nis, score_name='nis_score')\n",
    "sc.tl.score_genes(adata, nonReprog1, score_name='nonReprog1_score')\n",
    "sc.tl.score_genes(adata, nonReprog2, score_name='nonReprog2_score')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.draw_graph(adata, color=['cell_type', 'leiden_r0.75','mix_score'], palette='Set2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.draw_graph(adata, color=['cell_type', 'leiden_r0.75','early_prime_score'], palette='Set2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.draw_graph(adata, color=['cell_type', 'leiden_r0.75','primed_score'], palette='Set2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.draw_graph(adata, color=['cell_type', 'leiden_r0.75','naive_score'], palette='Set2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.draw_graph(adata, color=['cell_type', 'leiden_r0.75','nis_score'], palette='Set2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.draw_graph(adata, color=['cell_type', 'leiden_r0.75','nonReprog1_score'], palette='Set2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.draw_graph(adata, color=['cell_type', 'leiden_r0.75','nonReprog2_score'], palette='Set2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Lets look at everything together\n",
    "sc.pl.draw_graph(adata, color=['ANPEP','NANOG','ZIC2','DNMT3L','GATA3','fibro_score',\n",
    "                               'early_prime_score','primed_score','nis_score','naive_score','cell_type','leiden_r0.75'], palette='Set2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how our data looks compared to what is published"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "Image(\"ExtData_fig9e.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point we could run the marker gene scoring or take the top \"x\" genes for each cluster and do an enrichment analysis. However, I think that we are ready to name the clusters using the information that we have."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_cluster_names = ['early_primed', 'naive1', 'naive2', 'naive3', 'fibro','nonReprog2','fibro_nonReprog1',\n",
    "                     'primed','fibro_inter','Novel_Inter','mix','nonReprog1','12','naive4','14']\n",
    "adata.rename_categories('leiden_r0.75', new_cluster_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to export the image we would do that like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.tsne(adata, color=['leiden_r0.75'], legend_loc='on data' ,title='', frameon=False, \n",
    "           legend_fontsize=10, alpha=0.8, size=20, save = 'tsne_renamed_res0_75_winter2022.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.3 Subclustering "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To build on the basis clustering, we can now subcluster parts of the data to identify substructure within the identified cell types. Here, we subcluster the 'Novel Intermediate' population to see if we can find any difference between these cells in this cluster.\n",
    "\n",
    "Subclustering is normally performed at a lower resolution than on the entire dataset given that clustering is more sensitive when performed on a small subset of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Subcluster novel intermediates\n",
    "sc.tl.leiden(adata, restrict_to=('leiden_r0.75', ['Novel_Inter']), resolution=0.15, key_added='leiden_r0.75_NIS_sub')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Show the new clustering\n",
    "if 'leiden_r0.75_NIS_sub_colors' in adata.uns:\n",
    "    del adata.uns['leiden_r0.75_NIS_sub_colors']\n",
    "\n",
    "sc.pl.umap(adata, color='leiden_r0.75_NIS_sub')\n",
    "sc.pl.umap(adata, color='nis_score')\n",
    "sc.pl.umap(adata, color='cell_type')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The novel intermediate cluster broke up into three clusters and now we can find the marker genes in these clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Get the new marker genes\n",
    "sc.tl.rank_genes_groups(adata, groupby='leiden_r0.75_NIS_sub', key_added='rank_genes_r0.75_NIS_sub')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot the new marker genes\n",
    "sc.pl.rank_genes_groups(adata, key='rank_genes_r0.75_NIS_sub', groups=['Novel_Inter,0','Novel_Inter,1','Novel_Inter,2'], fontsize=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cluster 2 looks like it has some of the same markers seen in a similar report Castel et al. Cell Reports 2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "Image(\"Castel_et_al_2020_figS5.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obs['leiden_r0.75_NIS_sub'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.4 Trajectory Inference and pseudotime analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As our data set contains differentiation processes, we can investigate the differentiation trajectories in the data. This analysis is centred around the concept of 'pseudotime'. In pseudotime analysis a latent variable is inferred based on which the cells are ordered. This latent variable is supposed to measure the differentiation state along a trajectory.\n",
    "\n",
    "Pseudotime analysis is complicated when there are multiple trajectories in the data. In this case, the trajectory structure in the data must first be found before pseudotime can be inferred along each trajectory. The analysis is then called 'trajectory inference'.\n",
    "\n",
    "Once the pseudotime variable is inferred, we can test for genes that vary continuously along pseudotime. These genes are seen as being associated with the trajectory, and may play a regulatory role in the potential differentiation trajectory that the analysis found.\n",
    "\n",
    "Here, we measure the trajectory of the B/T/NK cells, which are supposed to be derived from a common progenitor. We also investigate which genes vary along pseudotime.\n",
    "\n",
    "Based on a recent comparison of pseudotime methods [Saelens et al., 2018], we have selected two of the top performing 'Monocle2', and 'Diffusion Pseudotime (DPT)'. Two methods were chosen as trajectory inference is a complex problem which is not yet solved. Different methods perform well on different types of trajectories. For example, 'Slingshot' was the top performer for simple bifurcating and multifurcating trajectories but takes more time and therefore we will skip this time; 'Monocle2' performed best for complex tree structures, and 'DPT' performed well in bifurcating trajectories. As the complexity of trajectories are generally not known, it is adviseable to compare trajectory inference outputs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Often we would first subset the data to include only the clusters that we care about. After subsetting it is important to recalculate the dimensionality reduction methods such as PCA, and diffusion maps, as the variability of the subsetted data set will be projected onto different basis vectors.\n",
    "\n",
    "Note that we subset the data to include only the cells we want to examine. Trajectory inference, and especially measuring gene expression changes over pseudotime can be a computationally expensive process, thus we often work with reduced gene sets that are informative of the variance in the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#first I am going to make a copy\n",
    "adata_copy = adata.copy()\n",
    "#Subset to highly variable genes\n",
    "sc.pp.highly_variable_genes(adata_copy, flavor='cell_ranger', n_top_genes=2500, subset=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we have subsetted the data to include only cell types that we assume are of interest, we recalculate the dimension reduction algorithms on this data. This ensures that for example the first few PCs capture only the variance in this data and not variance in parts of the full data set we have filtered out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Recalculating PCA for subset\n",
    "sc.pp.pca(adata_copy, svd_solver='arpack')\n",
    "sc.pl.pca(adata_copy)\n",
    "sc.pl.pca_variance_ratio(adata_copy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trajectory inference is often performed on PCA-reduced data, as is the case for Slingshot and Monocle2. To assess how many principal components (PCs) should be included in the low-dimensional representation we can use the 'elbow method'. This method involves looking for the 'elbow' in the plot of the variance ratio explained per PC. Above we can see the elbow at PC6. Thus the first six PCs are included in the slingshot data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_copy.obsm['X_pca'] = adata_copy.obsm['X_pca'][:,0:6]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.4.1 Monocle2 -- Skip This!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Monocle2 is written in R. The integration of R in this notebook is again achieved via the rpy2 interface. We use a specifically developed package called anndata2ri (https://www.github.com/flying-sheep/anndata2ri), that takes care of the conversion from an AnnData object to R. It should be noted that the convention for scRNA-seq data matrix storage in R is opposite to python. In R the expression matrix is stored as genes x cells rather than cells x genes. Thus, the matrix must be transposed before being input into the R function. This is already taken care of by anndata2ri.\n",
    "\n",
    "The monocle toolbox is one of the most popular trajectory inference toolboxes in the field. It has been shown to perform particularly well on complex trajectories. Although we only have a simple bi- or multifurcating trajectory here we use Monocle2 to confirm the suspected trajectory.\n",
    "\n",
    "Monocle2 is not as much a standalone trajectory inference method, than a trajectory inference toolbox. This toolbox includes preprocessing steps and a dimensionality reduction function. It is designed in such a way that integrating preprocessed data into the data structure is non-trivial. This is likely by intent to not interfere with the optimized pipeline. Thus, we use Monocle2 on our data set from the raw data. Our cluster labels are later overlayed onto the Monocle2 output to compare the trajectory with expectations.\n",
    "\n",
    "Implementation note:\n",
    "- this section closely follows the monocle 2 tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we need to import R for using in our python notebook and import all of the R packages that we need"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import rpy2.rinterface_lib.callbacks\n",
    "# import logging\n",
    "\n",
    "# from rpy2.robjects import pandas2ri\n",
    "# import anndata2ri\n",
    "# import rpy2.rinterface_lib.callbacks\n",
    "# import logging\n",
    "\n",
    "# from rpy2.robjects import pandas2ri\n",
    "# import anndata2ri\n",
    "# # Automatically convert rpy2 outputs to pandas dataframes\n",
    "# pandas2ri.activate()\n",
    "# anndata2ri.activate()\n",
    "# %load_ext rpy2.ipython\n",
    "\n",
    "# plt.rcParams['figure.figsize']=(8,8) #rescale figures\n",
    "# #sc.set_figure_params(dpi=200, dpi_save=300)\n",
    "# sc.logging.print_header()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%R\n",
    "\n",
    "# # Load all the R libraries we will be using in the notebook\n",
    "# library(RColorBrewer)\n",
    "# library(monocle) \n",
    "# library(ggplot2)\n",
    "# library(plyr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #Preprocessing for monocle\n",
    "# data_mat_mon = adata_copy.layers['counts'].T\n",
    "# var_mon=adata_copy.var.copy()\n",
    "# obs_mon=adata_copy.obs.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%R -i data_mat_mon -i obs_mon -i var_mon\n",
    "\n",
    "# #Set up the CellDataSet data structure\n",
    "# pd <- AnnotatedDataFrame(data = obs_mon)\n",
    "# fd <- AnnotatedDataFrame(data = var_mon)\n",
    "# colnames(data_mat_mon) <- rownames(pd)\n",
    "# rownames(data_mat_mon) <- rownames(fd)\n",
    "# ie_regions_cds <- newCellDataSet(cellData=data_mat_mon, phenoData=pd, featureData=fd, expressionFamily=negbinomial.size())\n",
    "\n",
    "# #Normalize the count data\n",
    "# ie_regions_cds <- estimateSizeFactors(ie_regions_cds)\n",
    "\n",
    "# #Calculate dispersions to filter for highly variable genes\n",
    "# ie_regions_cds <- estimateDispersions(ie_regions_cds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the plots above, it looks like most of the novel intermediates come from the TSC cell type, but what we have now labeled as \"naive\". Therefore, we are going to filter for the naive cells, novel inter, and mix clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%R\n",
    "\n",
    "# #Filter for naive, mix, and novel intermediates\n",
    "# cell_types = as.character(pData(ie_regions_cds)$leiden_r0.75)\n",
    "# cell_mask = rep(FALSE, length(cell_types))\n",
    "# cells_to_keep = c(\"naive1\", \"naive2\", \"naive3\", \"mix\", \"Novel_Inter\")\n",
    "# for (item in cells_to_keep) {cell_mask = cell_mask | startsWith(cell_types, item)}\n",
    "# print(\"Number of cells after filtering:\")\n",
    "# print(sum(cell_mask))\n",
    "# print(\"\")\n",
    "\n",
    "# #Filter highly variable genes from our analysis\n",
    "# hvg_mask = fData(ie_regions_cds)$highly_variable\n",
    "# ie_regions_cds <- ie_regions_cds[hvg_mask, cell_mask]\n",
    "\n",
    "# #Do dimensionality reduction\n",
    "# ie_regions_cds <- reduceDimension(ie_regions_cds, norm_method = 'vstExprs', reduction_method='DDRTree', verbose = F, max_components = 7)\n",
    "\n",
    "# #Run for the first time to get the ordering\n",
    "# ie_regions_cds <- orderCells(ie_regions_cds)\n",
    "\n",
    "# #Find the correct root state that corresponds to the 'naive1' cluster potentially\n",
    "# tab1 <- table(pData(ie_regions_cds)$State, pData(ie_regions_cds)$leiden_r0.75)\n",
    "# id = which(colnames(tab1) == 'naive1')\n",
    "# root_name = names(which.max(tab1[,id]))\n",
    "\n",
    "# #Run a second time to get the correct root state that overlaps with naive1 cells\n",
    "# ie_regions_cds <- orderCells(ie_regions_cds, root_state=root_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%R -w 1000 -h 800\n",
    "\n",
    "# #Get a nice colour map\n",
    "# custom_colour_map = brewer.pal(length(unique(pData(ie_regions_cds)$leiden_r0.75)),'Paired')\n",
    "\n",
    "# #Find the correct root state that coresponds to the 'naive1' cluster\n",
    "# tab1 <- table(pData(ie_regions_cds)$State, pData(ie_regions_cds)$leiden_r0.75)\n",
    "# id = which(colnames(tab1) == 'naive1')\n",
    "# root_name = names(which.max(tab1[,id]))\n",
    "\n",
    "# # Visualize with our cluster labels\n",
    "# options(repr.plot.width=5, repr.plot.height=4)\n",
    "# plot_complex_cell_trajectory(ie_regions_cds, color_by = 'leiden_r0.75', show_branch_points = T, \n",
    "#                              cell_size = 2, cell_link_size = 1, root_states = c(root_name)) +\n",
    "# scale_size(range = c(0.2, 0.2)) +\n",
    "# theme(legend.position=\"left\", legend.title=element_blank(), legend.text=element_text(size=rel(1.5))) +\n",
    "# guides(colour = guide_legend(override.aes = list(size=6))) + \n",
    "# scale_color_manual(values = custom_colour_map)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%R -w 600 -h 800\n",
    "\n",
    "# #Visualize pseudotime found\n",
    "# options(repr.plot.width=5, repr.plot.height=4)\n",
    "# plot_cell_trajectory(ie_regions_cds, color_by=\"Pseudotime\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.4.2 Diffusion Pseudotime (DPT) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We include Diffusion Pseudotime in the analysis to further support the found trajectories. Diffusion pseudotime is integrated into scanpy and is therefore easy to use with the current setup.\n",
    "\n",
    "DPT is based on diffusion maps, thus a diffusion map representation must be calculated prior to pseudotime inference. This in turn is based on a KNN graph embedding obtained from sc.pp.neighbors()."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.neighbors(adata_copy)\n",
    "sc.tl.diffmap(adata_copy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.diffmap(adata_copy, components='1,2', color='leiden_r0.75')\n",
    "sc.pl.diffmap(adata_copy, components='1,3', color='leiden_r0.75')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In DPT we must assign a root cell to infer pseudotime. In the plots we can observe that the most appropriate root will be either the naive1 or the nonReprog1 cells. We will try the nonReprog1 cell because it has the minimum DC1 and the maximum DC2 value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Find the CD4 T cell with the highest DC3 value to act as root for the diffusion pseudotime and compute DPT\n",
    "stem_mask = np.isin(adata_copy.obs['leiden_r0.75'], 'nonReprog1')\n",
    "max_stem_id = np.argmin(adata_copy.obsm['X_diffmap'][stem_mask,3])\n",
    "root_id = np.arange(len(stem_mask))[stem_mask][max_stem_id]\n",
    "adata_copy.uns['iroot'] = root_id\n",
    "\n",
    "#Compute dpt\n",
    "sc.tl.dpt(adata_copy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Visualize pseudotime over differentiation\n",
    "sc.pl.diffmap(adata_copy, components='1,3', color='dpt_pseudotime')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3.5 Partition-based graph abstraction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Partition-based graph abstraction (PAGA) is a method to reconcile clustering and pseudotemporal ordering. It can be applied to an entire dataset and does not assume that there are continuous trajectories connecting all clusters.\n",
    "\n",
    "As PAGA is integrated into scanpy, we can easily run it on the entire data set. Here we run and visualize PAGA with different clustering inputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.paga(adata, groups='leiden_r0.75')\n",
    "sc.pl.paga_compare(adata,fontsize=8)\n",
    "sc.pl.paga(adata,fontsize=8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The close connectivity of the two different subclusters shows the similarity of these cells. This result shows that not all connections in PAGA represent differentiation trajectories, but instead transcriptional similarity between states. Thus, further experiments are required to confirm potential lineage trajectories obtained via PAGA or other trajectory inference methods.\n",
    "\n",
    "We can do the same visualization on a umap layout."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.paga_compare(adata, basis='umap', fontsize=8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig1, ax1 = plt.subplots()\n",
    "sc.pl.umap(adata, size=40, ax=ax1, show=False)\n",
    "sc.pl.paga(adata, pos=adata.uns['paga']['pos'], show=False, node_size_scale=5, node_size_power=1, ax=ax1, text_kwds={'alpha':0})\n",
    "#plt.savefig('./figures/umap_paga_overlay_PBMCs.pdf', dpi=300, format='pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Implementation note:\n",
    "\n",
    "Note that the above plotting function only works when sc.pl.paga_compare(adata, basis='umap') is run before. The sc.pl.paga_compare() function stores the correct positions in adata.uns['paga']['pos'] to overlay the PAGA plot with a umap representation. To overlap PAGA with other representation, you can run sc.pl.paga_compare() with other basis parameters before plotting the combined plot."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regressing out the cell cycle effect will likely change how 'TA' cells are included in the trajectory. In this manner trajectory inference and graph abstraction can be iteratively improved. Trajectory inference and PAGA can be iteratively improved to better represent the biology. Knowing when to stop attempting to improve, or assessing when all of the relevant technical covariates have been taken into account can only be achieved with sufficient knowledge of the biological system, experience, and possibly some luck.\n",
    "\n",
    "It should also be noted that while an abstracted graph or an inferred trajectory can help to infer a lineage tree, experimental validation is necessary. Key driving forces in lineage specification might be lowly expressed genes and therefore neglected in the graph or even excluded in the HVG filtering."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.6 Write to file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Write to file\n",
    "adata.write('../Liu_et_al_example_winter2022.h5ad')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4.0 Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case study we went through the typical steps of an scRNA-seq data analysis workflow. We started with general preprocessing steps, which included cell and gene quality control, normalization, batch correction, selection of highly variable genes, visualization, and cell cycle scoring. In these steps the overall structure of the data is explored and filtered to produce optimal downstream analysis results. In the downstream scRNA-seq analysis section we then used methods to interpret the data and investigate particular parts of it. These steps included clustering and cluster identification via marker genes, trajectory inference, inferring an abstracted graph to relate clusters and trajectory inference in one visualization, and an example of differential expression although we did not have any experimental conditions to test.\n",
    "\n",
    "Any individual data analysis script will not always follow all of these steps. For example, trajectory inference may not be relevant where no differentiation processes are captured (this is probably the case with our data). Furthermore, a typical analysis will not always traverse the above steps as linearly as shown here. We have attempted to show this by suggesting how one can move back to data preprocessing to improve trajectory inference or PAGA. Indeed, all downstream analyses may require going back to tweak or add preprocessing steps to improve the downstream analysis results.\n",
    "\n",
    "A further aspect we stressed in this case study is that one should take care of which stage of data is being used to perform different analysis steps. Uncorrected, normalized data (or raw counts depending on the tool) are generally used for statistical tests where technical covariates such as batch can be included; corrected data should be used to visually compare results as the human eye will not be able to take into account technical differences that confound the comparison. When computing tests such as differential expression or finding marker genes, it is adviseable to use the full data set that is not restricted to find highly variable genes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('\\n'.join(f'{m.__name__} {m.__version__}' for m in globals().values() if getattr(m, '__version__', None)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aside: Using SCRAN on your own data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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. In fact, we don't estimate the exact number of molecules, but instead estimate cell-specific factors that should be proportional to the true number of molecules. These are called size factors. Normalized expression values are calculated by dividing the measured counts by the size factor for the cell.\n",
    "\n",
    "Based on comparisons of normalization methods for differential testing, and batch correction [Buettner et al., 2019], the normalization method implemented in the scran package performs best. This method requires a coarse clustering input to improve size factor esimation performance. Thus, we use a simple preprocessing approach and cluster the data at a low resolution to get an input for the size factor estimation. The basic preprocessing includes assuming all size factors are equal (library size normalization to counts per million - CPM) and log-transforming the count data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we need to import R for using in our python notebook and import all of the R packages that we need"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "scanpy_2021_test",
   "language": "python",
   "name": "scanpy_2021_test"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}