scATACseq_Tutorial4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Single-cell ATAC-seq Module 4\n",
    "\n",
    "## Background\n",
    "This tutorial demonstrates the principle behind the analysis of spatially resolvable single-cell ATACseq data. The transposase-accessible chromatin using sequencing (ATAC-seq) assay has been widely adopted since its development by Jason Buenrostro and colleagues in the Greenleaf lab in 2013, and it’s now one of the most popular approaches to investigate chromatin openness and its dynamics.\n",
    "\n",
    "ATAC-seq is used to identify regions of the genome that have open chromatin states that are generally associated with sites undergoing active transcription. Transcription factor binding sites and positions of nucleosomes can also be identified from the analysis of ATAC-seq data, potentially allowing important genetic pathways in the samples to be elucidated.\n",
    "\n",
    "The below graphic shows an overview of a typical ATAC sequencing experiment. These series of steps allow for discovery of chromatin accessibility dynamics. The flow displayed is: Batch correction -> Data transformation -> Dimension reduction -> Clustering -> Cell identity annotation.\n",
    "\n",
    "<img src=\"https://ars.els-cdn.com/content/image/1-s2.0-S2001037020303019-gr1.jpg\" alt=\"Drawing\" style=\"width:1000px;\"/>\n",
    "\n",
    "The above image is obtained from: Seungbyn Baek, Insuk Lee, Single-cell ATAC sequencing analysis: From data preprocessing to hypothesis generation, Computational and Structural Biotechnology Journal, Volume 18, 2020, Pages 1429-1439, ISSN 2001-0370,\n",
    "https://doi.org/10.1016/j.csbj.2020.06.012. PMID: [32637041](https://pubmed.ncbi.nlm.nih.gov/32637041/)\n",
    "\n",
    "## Standard single-cell ATAC-seq workflow\n",
    "The steps below encompass the standard preprocessing workflow. These represent the selection and filtration of cells based on QC metrics, data normalization and scaling, and the detection of highly variable features. \n",
    "\n",
    "The image below illustrates steps in a single-cell ATAC-Seq analysis workflow. Different from bulk tissue ATAC-Seq, single-cell ATAC-Seq can differentiate chromatin accessibility patterns in specific cell types. Firstly, reads are processed for quality control including steps such as demultiplexing, read alignment filtering, and duplicate removal. Normalization and dimensionality reduction are applied to address differences in sequencing depths and to get rid of the noise before clustering. After clusters are annotated, peak calling allows us to identify areas of the genome enriched for reads. This step is key to the overall goal of ATAC-Seq, identifying regions of accessible chromatin. Finally, we identify the differences in peaks among conditions in each cell cluster. Motif analysis allows for transcription factor discovery for the differentially enriched regions. \n",
    "\n",
    "<img src=\"https://media.springernature.com/full/springer-static/image/art%3A10.1038%2Fs41467-021-21583-9/MediaObjects/41467_2021_21583_Fig1_HTML.png?as=webp\" alt=\"Drawing\" style=\"width:1000px;\"/>\n",
    "\n",
    "The above image was obtained from the following paper: Fang, R., Preissl, S., Li, Y. et al. Comprehensive analysis of single-cell ATAC-seq data with SnapATAC. Nat Commun 12, 1337 (2021). https://doi.org/10.1038/s41467-021-21583-9.\n",
    "PMID [33637727](https://pubmed.ncbi.nlm.nih.gov/33637727/)\n",
    "\n",
    "This tutorial will cover the following tasks. Prior to these steps, sample preparation, library construction and standard NGS short-read sequencing for massively parallel epigenomic profiling of thousands of individual cells will be performed.\n",
    "- QC and selecting cells\n",
    "- Normalization\n",
    "- Identifying highly variable features\n",
    "- Scaling the data\n",
    "- Linear Dimensional reduction \n",
    "- Determining Dimensionality\n",
    "- Clustering\n",
    "- Assigning cell type identity to clusters\n",
    "- Detecting spatially variable features\n",
    "- Interactive visualization\n",
    "- Integration with single-cell RNA-seq data\n",
    "\n",
    "## Data and Analysis\n",
    "The dataset we use here comes from Lareau et al., Nat Biotech 2019(PMID: [33637727](https://pubmed.ncbi.nlm.nih.gov/33637727/)), one of the highest throughput single-cell ATAC-seq experiments to date. They assay 136,000 resting and stimulated bone marrow-derived cells and study the different cell types and lineages that emerge. The paper analyzes resting and stimulated cells separately. We also find that these two populations split into largely distinct clusters, so here we only focus on the 60,000 resting cells from this experiment.\n",
    "\n",
    "The cell metadata from this experiment includes annotated putative cell types based on chromVAR, Louvain clustering, and downstream cluster analysis. We compare our de novo identified clusters to these annotated clusters from the paper (and find good agreement).\n",
    "\n",
    "The peak-cell matrix along with peak and cell metadata files are publicly available on GEO, [GSE123580](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE123580). We processed the dataset to include only cells in the 'Resting' condition and peaks with nonzero coverage. In this tutorial, you will download (1) the processed peak-cell count matrix for this dataset (.h5ad), (2) the set of nonzero peak names (.npy), and (3) the cell metadata (.csv) from a bucket. In addition, there are 2 custom function files that need to be imported for the notebook to execute properly: 1) [rapids_scanpy_funcs.py](Tutorial4/rapids_scanpy_funcs.py) and 2) [utils.py](Tutorial4/utils.py). \n",
    "\n",
    "#### GPU Acceleration\n",
    "\n",
    "Single-cell analysis is complex and compute intensive leading to long analysis times on an average HPC. Using GPUs speeds up the analysis. Google Cloud provides images with the ability for NVIDIA GPU drivers to be auto installed and some pre-requisite software pre-installed. The details are highlighted in the [README](README.md) of this module. \n",
    "\n",
    "This tutorial was adapted from the following NVIDIA GitHub repository (https://github.com/NVIDIA-Genomics-Research/rapids-single-cell-examples/blob/master/README.md). Follow the provided GitHub link to view acceleration and cost statistics from runs by the NVIDIA research team from December of 2020. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Take quiz - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install jupyterquiz==2.0.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install jupytercards"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import jupyterquiz\n",
    "from jupyterquiz import display_quiz\n",
    "display_quiz(\"quiz_files/Quiz-1.json\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert-info\" style=\"font-size:200%\">\n",
    "STEP 1: Set Up Environment\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Installing gcloud, including gsutil and other google command line tools\n",
    "!chmod +x Tutorial4/install-gcloud.sh\n",
    "!bash Tutorial4/install-gcloud.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import Sample Files\n",
    "Now lets create some folders to stay organized and copy over our prepared raw files. We're going to create a directory called \"Tutorial4\" which we'll use for this module. We'll then create sub-folders for our input files and for the files that we'll be creating during this module."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following files will be downloaded when you execute the below command\n",
    "- processed peak-cell count matrix for this dataset (.h5ad)\n",
    "- set of nonzero peak names (.npy)\n",
    "- cell metadata (.csv)\n",
    "are stored inside the data folder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# These commands create our directory structure.\n",
    "#!cd $HOMEDIR\n",
    "!mkdir -p Tutorial4/InputFiles\n",
    "#!echo $PWD\n",
    "\n",
    "# These commands help identify the Google Cloud storage bucket where the example files are held.\n",
    "original_bucket = \"gs://nigms-sandbox/unmc_atac_data_examples/Tutorial4\"\n",
    "\n",
    "# This command copies our example files to the Tutorial4/InputFiles folder that we created above.\n",
    "!gsutil cp $original_bucket/dsci_resting_cell_metadata.csv Tutorial4/InputFiles\n",
    "!gsutil cp $original_bucket/dsci_resting_nonzeropeaks.h5ad Tutorial4/InputFiles\n",
    "!gsutil cp $original_bucket/dsci_resting_peaknames_nonzero.npy Tutorial4/InputFiles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Hierarchical Data Format version 5 (H5AD/HDF5), is an open source file format that supports large, complex, heterogeneous data. Cell Ranger (a set of analysis pipelines) generates this output file with per-molecule information in HDF5 format."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's make sure that the files copied correctly. \n",
    "### You should see three files after running the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls Tutorial4/InputFiles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Installing required packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#commented packages are already in the image to see full list see pip_packages_preinstalled.txt\n",
    "#!pip install numpy\n",
    "!pip install scanpy\n",
    "!pip install anndata\n",
    "#!pip install matplotlib\n",
    "#!pip install rmm\n",
    "#!pip install cupy \n",
    "#!pip install cudf\n",
    "!pip install utils\n",
    "#!pip install rapids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# *Restart Kernel here*\n",
    "Do the following -> Menu tab > Kernel > Restart Kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#did not restart kernel or do this step as its already done\n",
    "!pip install scanpy wget"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# *Restart Kernel here again*\n",
    "Do the following -> Menu tab > Kernel > Restart Kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#did not run above restart kernel\n",
    "#got rid of , wget after import os\n",
    "import numpy as np\n",
    "import scanpy as sc\n",
    "import anndata\n",
    "import matplotlib.pyplot as plt\n",
    "from collections import OrderedDict\n",
    "import time\n",
    "import os\n",
    "import rmm\n",
    "import cudf\n",
    "import cupy as cp\n",
    "from cuml.decomposition import PCA\n",
    "from cuml.manifold import TSNE\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore', 'Expected ')\n",
    "warnings.simplefilter('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#installing custom python functions by nvidia research team \n",
    "import Tutorial4.rapids_scanpy_funcs as rapids_scanpy_funcs\n",
    "import Tutorial4.utils as utils\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmm.reinitialize(\n",
    "    managed_memory=True, # Allows oversubscription\n",
    "    devices=0, # GPU device IDs to register. By default registers only GPU 0.\n",
    ")\n",
    "\n",
    "cp.cuda.set_allocator(rmm.rmm_cupy_allocator)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert-info\" style=\"font-size:200%\">\n",
    "STEP 2: Setting Up Parameters\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### In the below set of parameters, you will begin defining each constant to enable successful calling of peaks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filtering peaks\n",
    "n_top_peaks = 5000 # Number of peaks to retain\n",
    "\n",
    "# PCA\n",
    "n_components = 50 # Number of principal components to compute\n",
    "\n",
    "# t-SNE\n",
    "tsne_n_pcs = 20 # Number of principal components to use for t-SNE\n",
    "\n",
    "# KNN\n",
    "n_neighbors = 10 # Number of nearest neighbors for KNN graph\n",
    "knn_n_pcs = 50 # Number of principal components to use for finding nearest neighbors\n",
    "\n",
    "# UMAP\n",
    "umap_min_dist = 0.4 \n",
    "umap_spread = 1.0\n",
    "\n",
    "# Differential peaks\n",
    "n_diff_peaks = 50 # Number of differential peaks to compute for each cluster\n",
    "\n",
    "# Marker genes\n",
    "# Gene coordinates are in GRCh37, from grch37.ensembl.org\n",
    "markers = ['MS4A1', 'LEF1', 'NKG7', 'TREM1', 'GATA1', 'PAX5']\n",
    "marker_coords = OrderedDict()\n",
    "marker_coords['MS4A1'] = ['chr11', 60223225, 60238233, 'fwd']\n",
    "marker_coords['LEF1'] = ['chr4', 108968701, 109090112, 'rev']\n",
    "marker_coords['NKG7'] = ['chr19', 51874860, 51875969, 'rev']\n",
    "marker_coords['TREM1'] = ['chr6', 41235664, 41254457, 'rev']\n",
    "marker_coords['GATA1'] = ['chrX', 48644962, 48652716, 'fwd']\n",
    "marker_coords['PAX5'] = ['chr9', 36833272, 37034103, 'rev']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_time = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Input data\n",
    "\n",
    "In the cell below, we provide the path to the sparse .h5ad file containing the peak x cell count matrix to analyze."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "input_file = \"Tutorial4/InputFiles/dsci_resting_nonzeropeaks.h5ad\"\n",
    "\n",
    "adata = anndata.read_h5ad(input_file)\n",
    "adata_raw = adata.copy()\n",
    "print(adata.X.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flash Card - 1: Learn the primary advantages of single cell ATAC-seq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import jupytercards\n",
    "from jupytercards import display_flashcards\n",
    "display_flashcards(\"quiz_files/my_cards-1.json\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert-info\" style=\"font-size:200%\">\n",
    "STEP 3: Preprocessing\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "preprocessing_start = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TF-IDF normalization\n",
    "\n",
    "We use log1p(pseudocount / tf)*idf. It uses the raw binary count matrix as the TF matrix rather than dividing by the total reads per cell, which obviates the need for log scaling. The resulting TF-IDF transformed matrix will be used as input into PCA without any further normalization or scaling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "lognormalized = utils.logtf_idf(adata.X, pseudocount=10**3)\n",
    "adata.X = lognormalized"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Frequency-based peak selection\n",
    "\n",
    "We want to choose peaks that have highest discriminative power between cells. A good proxy for discriminative power is variance; since peak data (especially on a sparse dataset like this) is essentially binary, the highest variance peaks are the ones with closest to 50% representation across cells (p=1/2 maximizes Bernoulli variance). The most frequent peak in this dataset only has 33% representation, so in practice we just end up taking the top N most frequent peaks, where N is a parameter we set above (default 5000)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "adata = utils.filter_peaks(adata, n_top_peaks)\n",
    "print(adata.X.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Preprocessing time: %.2fsec\" % (time.time() - preprocessing_start))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert-info\" style=\"font-size:200%\">\n",
    "STEP 4: Clustering & Visualization\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clustering is the task of dividing the cells into a number of groups such that data points (open chromatin peaks) in the same groups are more similar to other data points in the same group and dissimilar to the data points in other groups."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Principal component analysis, or PCA, is a dimensionality-reduction method that is often used to reduce the dimensionality of large data sets, by transforming a large set of variables into a smaller one that still contains most of the information in the large set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "adata = anndata.AnnData(X=adata.X.todense(),\n",
    "                        obs=adata.obs,\n",
    "                        var=adata.var)\n",
    "adata.obsm[\"X_pca\"] = PCA(n_components=n_components).fit_transform(adata.X)\n",
    "adata.obsm[\"X_pca\"].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flash Card - 2: Learn the benefits of single-cell ATAC-seq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import jupytercards\n",
    "from jupytercards import display_flashcards\n",
    "display_flashcards(\"quiz_files/my_cards-2.json\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# UMAP visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uniform Manifold Approximation and Projection (UMAP) is a dimensionality reduction technique that can be used for visualization similarly to t-SNE, but also for general non-linear dimension reduction.\n",
    "\n",
    "UMAP algorithm is founded on three assumptions about the data:\n",
    "\n",
    "1. The data is uniformly distributed on a Riemannian manifold;\n",
    "2. The Riemannian metric is locally constant (or can be approximated as such);\n",
    "3. The manifold is locally connected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "sc.pp.neighbors(adata, n_neighbors=n_neighbors, n_pcs=knn_n_pcs, method='rapids')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "sc.tl.umap(adata, min_dist=umap_min_dist, spread=umap_spread, method='rapids')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flash Card - 3: Understanding the benefits of single cell ATAC-seq compared to bulk ATAC-seq?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import jupytercards\n",
    "from jupytercards import display_flashcards\n",
    "display_flashcards(\"quiz_files/my_cards-3.json\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Graph clustering\n",
    "Below, we show how to use the RAPIDS Louvain clustering or Leiden clustering functions to cluster the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "sc.tl.louvain(adata, flavor='rapids')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "adata.obs['leiden'] = rapids_scanpy_funcs.leiden(adata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we plot the UMAP visualization of the cells colored by the Louvain clusters. Next, we color the cells using the cell type annotations reported in the source paper. We see that the clusters we discovered match closely with the reported cell types."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "fig, axs = plt.subplots(nrows = 1, ncols = 3, figsize = (12, 3), dpi=100)\n",
    "\n",
    "sc.pl.umap(adata, color=['louvain'], ax=axs[0], show=False)\n",
    "sc.pl.umap(adata, color=['leiden'], ax=axs[1], show=False)\n",
    "sc.pl.umap(adata, color=['cell_type'], ax=axs[2], show=False)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Take quiz - 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display_quiz(\"quiz_files/Quiz-2.json\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## t-SNE visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "T-distributed neighbor embedding (t-SNE) is a dimensionality reduction technique that helps visualize high-dimensional data sets. It takes the original data and matches both distributions to determine how to best represent this data using fewer dimensions.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "adata.obsm['X_tsne'] = TSNE().fit_transform(adata.obsm[\"X_pca\"][:, :tsne_n_pcs])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "fig, axs = plt.subplots(nrows = 1, ncols = 3, figsize = (12, 3), dpi=100)\n",
    "\n",
    "sc.pl.tsne(adata, color=['louvain'], ax=axs[0], show=False)\n",
    "sc.pl.tsne(adata, color=['leiden'], ax=axs[1], show=False)\n",
    "sc.pl.tsne(adata, color=['cell_type'], ax=axs[2], show=False)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Marker genes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now compute a gene activity score representing the activity of each marker gene in each cell. We use the adata_raw object to compute this score. This allows us to consider all peaks, not just the subset we used for clustering and UMAP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# find which peaks overlap with each marker gene (+ 5kb upstream)\n",
    "peaks_per_gene = {gene : [] for gene in markers}\n",
    "\n",
    "for (i, peak) in enumerate(adata_raw.var_names):\n",
    "    chrom, coords = peak.strip().split(':')\n",
    "    start, end = [int(c) for c in coords.split('-')]\n",
    "    for gene in markers:\n",
    "        if utils.overlap(marker_coords[gene], [chrom, start, end], 5000, 0):\n",
    "            peaks_per_gene[gene].append((i, peak))\n",
    "\n",
    "for gene in peaks_per_gene:\n",
    "    print(f\"{gene} has {len(peaks_per_gene[gene])} nearby peaks\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Calculate per-cell gene activity scores for each marker gene\n",
    "cell_gene_scores = np.zeros((adata.X.shape[0], len(marker_coords.keys())))\n",
    "for (i, gene) in enumerate(peaks_per_gene.keys()):\n",
    "    indices = [x[0] for x in peaks_per_gene[gene]]\n",
    "    peak_scores = adata_raw.X[:, indices].toarray()\n",
    "    cell_gene_scores[:, i] = np.sum(peak_scores, axis=1) + 1\n",
    "\n",
    "# Store per-cell gene activity scores with adata\n",
    "for (i, gene) in enumerate(markers):\n",
    "    adata.obs[gene] = np.log2(cell_gene_scores[:, i].tolist())\n",
    "    adata.obs.loc[np.isneginf(adata.obs[gene]), gene] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Take quiz - 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display_quiz(\"quiz_files/Quiz-3.json\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now visualize the UMAP plot colored by activity of marker genes. We observe that the activation of cell-type specific marker genes correlates well with the cell types identified by Louvain clustering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "fig, axs = plt.subplots(nrows=2, ncols=3, figsize = (8, 4), dpi=100)\n",
    "axs = axs.flatten()\n",
    "for (i, gene) in enumerate(markers):\n",
    "    sc.pl.umap(adata, color_map=\"Blues\", color=gene, ax = axs[i], show=False)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert-info\" style=\"font-size:200%\">\n",
    "STEP 5: Finding differential peaks\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Find Differential peaks\n",
    "We perform an accelerated logistic regression-based differential peak computation using RAPIDS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "cluster_labels = cudf.Series.from_categorical(adata.obs[\"louvain\"].cat)\n",
    "var_names = cudf.Series(adata.var_names)\n",
    "dense_gpu_array = cp.sparse.csr_matrix(cp.array(adata.X))\n",
    "\n",
    "scores, names, reference = rapids_scanpy_funcs.rank_genes_groups(\n",
    "    dense_gpu_array,\n",
    "    cluster_labels, \n",
    "    var_names, \n",
    "    penalty='none',\n",
    "    n_genes=n_diff_peaks, groups='all', reference='rest')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "fig, axs = plt.subplots(len(names[0])//3 + 1, 3, figsize=(11,11), dpi=100)\n",
    "axs = axs.flatten()\n",
    "\n",
    "for (i, peak) in enumerate(names[0]):\n",
    "    sc.pl.umap(adata, color=peak, \n",
    "           ax=axs[i], show=False,\n",
    "           vmax=10, vmin=-0.5,\n",
    "           cmap='Blues'\n",
    "          )\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Full time: %.2fsec\" % (time.time() - start_time))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Congratulations! You have successfully completed the single-cell tutorial\n",
    "You successfully performed peak-cell matrix, then perform peak selection, normalization, dimensionality reduction, clustering, and visualization on single cell ATAC-seq data. You also visualized regulatory activity at marker genes and computed differential peaks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\">\n",
    "    <b>&#128721; Caution:</b> Remember to shut down your VM after you are finished with your work in order to avoid incurring additional charges.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 4
}