scanpy_tutorial_HoloViz.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ddc32c04",
   "metadata": {},
   "source": [
    "# [HoloViz](#HoloViz) stuff is at the bottom"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "404b48be-6e11-4539-8524-4213162d909c",
   "metadata": {},
   "source": [
    "# Preprocessing and clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1990713-5968-4fc7-aeb6-5a842606fd8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Core scverse libraries\n",
    "import scanpy as sc\n",
    "import anndata as ad\n",
    "\n",
    "# Data retrieval\n",
    "import pooch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "161629d6-b1b6-4faa-8bca-01a841e0ec10",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.settings.set_figure_params(dpi=50, facecolor=\"white\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "3211b168-0575-43fa-945c-405ee9438098",
   "metadata": {},
   "source": [
    "The data used in this basic preprocessing and clustering tutorial was collected from bone marrow mononuclear cells of healthy human donors and was part of [openproblem's NeurIPS 2021 benchmarking dataset](https://openproblems.bio/competitions/neurips_2021/) {cite}`Luecken2021`. The samples used in this tutorial were measured using the 10X Multiome Gene Expression and Chromatin Accessability kit. \n",
    "\n",
    "\n",
    "We are reading in the count matrix into an [AnnData](https://anndata.readthedocs.io/en/latest/tutorials/notebooks/getting-started.html) object, which holds many slots for annotations and different representations of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0448e97",
   "metadata": {},
   "outputs": [],
   "source": [
    "EXAMPLE_DATA = pooch.create(\n",
    "    path=pooch.os_cache(\"scverse_tutorials\"),\n",
    "    base_url=\"doi:10.6084/m9.figshare.22716739.v1/\",\n",
    ")\n",
    "EXAMPLE_DATA.load_registry_from_doi()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9bdf7f10-e9b9-4929-a7e2-7ac373d769d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "samples = {\n",
    "    \"s1d1\": \"s1d1_filtered_feature_bc_matrix.h5\",\n",
    "    \"s1d3\": \"s1d3_filtered_feature_bc_matrix.h5\",\n",
    "}\n",
    "adatas = {}\n",
    "\n",
    "for sample_id, filename in samples.items():\n",
    "    path = EXAMPLE_DATA.fetch(filename)\n",
    "    sample_adata = sc.read_10x_h5(path)\n",
    "    sample_adata.var_names_make_unique()\n",
    "    adatas[sample_id] = sample_adata\n",
    "\n",
    "adata = ad.concat(adatas, label=\"sample\")\n",
    "adata.obs_names_make_unique()\n",
    "print(adata.obs[\"sample\"].value_counts())\n",
    "adata"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c673fe14-df28-4d53-abf6-3a852040e571",
   "metadata": {},
   "source": [
    "The data contains ~8,000 cells per sample and 36,601 measured genes. We'll now investigate these with a basic preprocessing and clustering workflow."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e306dde3-b587-4561-b5ec-e8a2ac5e2c13",
   "metadata": {},
   "source": [
    "## Quality Control"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "cd32b625-7f24-4b9c-9564-dd97c57dc1b2",
   "metadata": {},
   "source": [
    "The scanpy function {func}`~scanpy.pp.calculate_qc_metrics` calculates common quality control (QC) metrics, which are largely based on `calculateQCMetrics` from scater {cite}`McCarthy2017`. One can pass specific gene population to {func}`~scanpy.pp.calculate_qc_metrics` in order to calculate proportions of counts for these populations. Mitochondrial, ribosomal and hemoglobin genes are defined by distinct prefixes as listed below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75df4b20-29ca-4e12-a878-94698a1b67d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# mitochondrial genes, \"MT-\" for human, \"Mt-\" for mouse\n",
    "adata.var[\"mt\"] = adata.var_names.str.startswith(\"MT-\")\n",
    "# ribosomal genes\n",
    "adata.var[\"ribo\"] = adata.var_names.str.startswith((\"RPS\", \"RPL\"))\n",
    "# hemoglobin genes\n",
    "adata.var[\"hb\"] = adata.var_names.str.contains(\"^HB[^(P)]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efd36ca1-75c7-4716-8766-a5571951537b",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.calculate_qc_metrics(\n",
    "    adata, qc_vars=[\"mt\", \"ribo\", \"hb\"], inplace=True, log1p=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b525671d-cfd3-4989-ad06-91c6edfa3af8",
   "metadata": {},
   "source": [
    "One can now inspect violin plots of some of the computed QC metrics:\n",
    "\n",
    "* the number of genes expressed in the count matrix\n",
    "* the total counts per cell\n",
    "* the percentage of counts in mitochondrial genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9db90e7-57d8-4f78-b1f9-ffc3698df091",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f8e3e58-a05a-4fd1-ab77-f3b14232dbf8",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.violin(\n",
    "    adata,\n",
    "    [\"n_genes_by_counts\", \"total_counts\", \"pct_counts_mt\"],\n",
    "    jitter=0.4,\n",
    "    multi_panel=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b78eec5-6b08-4d4f-9c54-a3688984b899",
   "metadata": {},
   "source": [
    "Additionally, it is useful to consider QC metrics jointly by inspecting a scatter plot colored by `pct_counts_mt`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ddfd195-10cd-419d-b984-39bc29dfe56b",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.scatter(adata, \"total_counts\", \"n_genes_by_counts\", color=\"pct_counts_mt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2de6f565-7a2a-4ff5-a783-9dab32a007a2",
   "metadata": {},
   "source": [
    "Based on the QC metric plots, one could now remove cells that have too many mitochondrial genes expressed or too many total counts by setting manual or automatic thresholds. However, sometimes what appears to be poor QC metrics can be driven by real biology so we suggest starting with a very permissive filtering strategy and revisiting it at a later point. We therefore now only filter cells with less than 100 genes expressed and genes that are detected in less than 3 cells. \n",
    "\n",
    "Additionally, it is important to note that for datasets with multiple batches, quality control should be performed for each sample individually as quality control thresholds can very substantially between batches. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61508aa7-31f1-403c-96c2-09460427bdae",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.filter_cells(adata, min_genes=100)\n",
    "sc.pp.filter_genes(adata, min_cells=3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "71c53c35-7f0d-4c8d-90cc-fdc7687e828d",
   "metadata": {},
   "source": [
    "### Doublet detection\n",
    "\n",
    "As a next step, we run a doublet detection algorithm. Identifying doublets is crucial as they can lead to misclassifications or distortions in downstream analysis steps. Scanpy contains the doublet detection method Scrublet {cite}`Wolock2019`. Scrublet predicts cell doublets using a nearest-neighbor classifier of observed transcriptomes and simulated doublets. {func}`scanpy.pp.scrublet` adds `doublet_score` and `predicted_doublet` to `.obs`. One can now either filter directly on `predicted_doublet` or use the `doublet_score` later during clustering to filter clusters with high doublet scores. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bc52329-32e3-4c04-a5a6-be4184391219",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.scrublet(adata, batch_key=\"sample\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85bcdb2e",
   "metadata": {},
   "source": [
    "We can remove doublets by either filtering out the cells called as doublets, or waiting until we've done a clustering pass and filtering out any clusters with high doublet scores."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "eca84d9a-1da7-40c5-ab50-a6bbc5cd6955",
   "metadata": {},
   "source": [
    ":::{seealso}\n",
    "Alternative methods for doublet detection within the scverse ecosystem are [DoubletDetection](https://github.com/JonathanShor/DoubletDetection) and [SOLO](https://docs.scvi-tools.org/en/stable/user_guide/models/solo.html). You can read more about these in the [Doublet Detection chapter](https://www.sc-best-practices.org/preprocessing_visualization/quality_control.html#doublet-detection) of Single Cell Best Practices.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31ae375a-5f9d-4ff6-a7b2-52d13d202fa8",
   "metadata": {},
   "source": [
    "## Normalization\n",
    "\n",
    "The next preprocessing step is normalization. A common approach is count depth scaling with subsequent log plus one (log1p) transformation. Count depth scaling normalizes the data to a “size factor” such as the median count depth in the dataset, ten thousand (CP10k) or one million (CPM, counts per million). The size factor for count depth scaling can be controlled via `target_sum` in `pp.normalize_total`. We are applying median count depth normalization with log1p transformation (AKA log1PF)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fe44a18",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Saving count data\n",
    "adata.layers[\"counts\"] = adata.X.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "053aed6d-f793-442d-afa8-46def6054b23",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normalizing to median total counts\n",
    "sc.pp.normalize_total(adata)\n",
    "# Logarithmize the data\n",
    "sc.pp.log1p(adata)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "1f634d44-4e23-4e0f-96c0-443b3a5ffda9",
   "metadata": {},
   "source": [
    "## Feature selection\n",
    "\n",
    "As a next step, we want to reduce the dimensionality of the dataset and only include the most informative genes. This step is commonly known as feature selection. The scanpy function `pp.highly_variable_genes` annotates highly variable genes by reproducing the implementations of Seurat {cite}`Satija2015`, Cell Ranger {cite}`Zheng2017`, and Seurat v3 {cite}`Stuart2019` depending on the chosen `flavor`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6241df3f-a23b-4f06-8c4e-a2de39200f05",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.highly_variable_genes(adata, n_top_genes=2000, batch_key=\"sample\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "654984e9-e511-4e32-b025-054d41a1aabd",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.highly_variable_genes(adata)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe4e74b2-c177-47d1-b659-041e8b796231",
   "metadata": {},
   "source": [
    "## Dimensionality Reduction\n",
    "Reduce the dimensionality of the data by running principal component analysis (PCA), which reveals the main axes of variation and denoises the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a4d36d8-3238-47a5-8950-049fd0473205",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.pca(adata)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "bb54e365-3989-432e-bce0-b0a8e46ac08c",
   "metadata": {},
   "source": [
    "Let us inspect the contribution of single PCs to the total variance in the data. This gives us information about how many PCs we should consider in order to compute the neighborhood relations of cells, e.g. used in the clustering function {func}`~scanpy.tl.leiden` or {func}`~scanpy.tl.tsne`. In our experience, there does not seem to be signifigant downside to overestimating the numer of principal components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b38f1a27-6c1b-4160-b28e-e172bb325d4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.pca_variance_ratio(adata, n_pcs=50, log=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa711f1f",
   "metadata": {},
   "source": [
    "You can also plot the principal components to see if there are any potentially undesired features (e.g. batch, QC metrics) driving signifigant variation in this dataset. In this case, there isn't anything too alarming, but it's a good idea to explore this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58f92649",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.pca(\n",
    "    adata,\n",
    "    color=[\"sample\", \"sample\", \"pct_counts_mt\", \"pct_counts_mt\"],\n",
    "    dimensions=[(0, 1), (2, 3), (0, 1), (2, 3)],\n",
    "    ncols=2,\n",
    "    size=2,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27147fca-c120-468d-8a5c-6479ece27848",
   "metadata": {},
   "source": [
    "## Nearest neighbor graph constuction and visualization\n",
    "\n",
    "Let us compute the neighborhood graph of cells using the PCA representation of the data matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce3adfa7-0be0-4232-88fa-e0540c5c3319",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pp.neighbors(adata)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a0f7719-e418-401a-b006-c1eb524d7a2e",
   "metadata": {},
   "source": [
    "This graph can then be embedded in two dimensions for visualiztion with UMAP (McInnes et al., 2018):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d0151a0-ce62-40e6-ba1b-2f4505b6b4b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.umap(adata)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "857f818f-66c7-40c7-b834-47cf5f02a1fa",
   "metadata": {},
   "source": [
    "We can now visualize the UMAP according to the `sample`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74aba38e-9c07-484b-a222-5f2ec108e7b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(\n",
    "    adata,\n",
    "    color=\"sample\",\n",
    "    # Setting a smaller point size to get prevent overlap\n",
    "    size=2,\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "6d472937-b5b8-4d38-88d5-2af4d0718f1e",
   "metadata": {},
   "source": [
    "Even though the data considered in this tutorial includes two different samples, we only observe a minor batch effect and we can continue with clustering and annotation of our data. \n",
    "\n",
    "If you inspect batch effects in your UMAP it can be beneficial to integrate across samples and perform batch correction/integration. We recommend checking out [`scanorama`](https://github.com/brianhie/scanorama) and [`scvi-tools`](https://scvi-tools.org) for batch integration."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "358da4af-c7ed-4479-87a7-ecd9d9765e1a",
   "metadata": {},
   "source": [
    "## Clustering\n",
    "\n",
    "As with Seurat and many other frameworks, we recommend the Leiden graph-clustering method (community detection based on optimizing modularity) {cite}`Traag2019`. Note that Leiden clustering directly clusters the neighborhood graph of cells, which we already computed in the previous section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11bb19be-3df1-4a17-ad84-16b50a9348ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using the igraph implementation and a fixed number of iterations can be significantly faster, especially for larger datasets\n",
    "sc.tl.leiden(adata, flavor=\"igraph\", n_iterations=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e616f3f6-69f8-4c66-ac00-48d328dcf755",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(adata, color=[\"leiden\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d7410a1-333a-4f21-b53e-91af768b07c7",
   "metadata": {},
   "source": [
    "## Re-assess quality control and cell filtering \n",
    "\n",
    "As indicated before, we will now re-assess our filtering strategy by visualizing different QC metrics using UMAP. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1818fbe-0e7d-467f-bdd1-3b690c32ad26",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(\n",
    "    adata,\n",
    "    color=[\"leiden\", \"predicted_doublet\", \"doublet_score\"],\n",
    "    # increase horizontal space between panels\n",
    "    wspace=0.5,\n",
    "    size=3,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91ad2b94-ab9b-48a3-aa0c-551203d3cf91",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(\n",
    "    adata,\n",
    "    color=[\"leiden\", \"log1p_total_counts\", \"pct_counts_mt\", \"log1p_n_genes_by_counts\"],\n",
    "    wspace=0.5,\n",
    "    ncols=2,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dece397a-e473-4633-b552-a7e430095c10",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Manual cell-type annotation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5e54a09",
   "metadata": {},
   "source": [
    ":::{note}\n",
    "This section of the tutorial is expanded upon using prior knowledge resources like automated assignment and gene enrichment in the scverse tutorial [here](https://scverse-tutorials.readthedocs.io/en/latest/notebooks/basic-scrna-tutorial.html#cell-type-annotation)\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ff5d831",
   "metadata": {},
   "source": [
    "Cell type annotation is laborous and repetitive task, one which typically requires multiple rounds of subclustering and re-annotation. It's difficult to show the entirety of the process in this tutorial, but we aim to show how the tools scanpy provides assist in this process."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc966c08-e87e-499a-9a9e-db40b243deb9",
   "metadata": {},
   "source": [
    "We have now reached a point where we have obtained a set of cells with decent quality, and we can proceed to their annotation to known cell types. Typically, this is done using genes that are exclusively expressed by a given cell type, or in other words these genes are the marker genes of the cell types, and are thus used to distinguish the heterogeneous groups of cells in our data. Previous efforts have collected and curated various marker genes into available resources, such as [CellMarker](http://bio-bigdata.hrbmu.edu.cn/CellMarker/), [TF-Marker](http://bio.liclab.net/TF-Marker/), and [PanglaoDB](https://panglaodb.se/). The [cellxgene gene expression tool](https://cellxgene.cziscience.com/gene-expression) can also be quite useful to see which cell types a gene has been expressed in across many existing datasets."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6bc8b4d-dd05-41c2-befd-048a5eea6afd",
   "metadata": {},
   "source": [
    "Commonly and classically, cell type annotation uses those marker genes subsequent to the grouping of the cells into clusters. So, let's generate a set of clustering solutions which we can then use to annotate our cell types. Here, we will use the Leiden clustering algorithm which will extract cell communities from our nearest neighbours graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "148a1050-a073-4d82-b568-3cb943f9fd34",
   "metadata": {},
   "outputs": [],
   "source": [
    "for res in [0.02, 0.5, 2.0]:\n",
    "    sc.tl.leiden(\n",
    "        adata, key_added=f\"leiden_res_{res:4.2f}\", resolution=res, flavor=\"igraph\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbffb2d3-e8fe-4bfb-aa8a-8d490254dcb1",
   "metadata": {},
   "source": [
    "Notably, the number of clusters that we define is largely arbitrary, and so is the `resolution` parameter that we use to control for it. As such, the number of clusters is ultimately bound to the stable and biologically-meaningful groups that we can ultimately distringuish, typically done by experts in the corresponding field or by using expert-curated prior knowledge in the form of markers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1275cb4-9f1e-4f37-8c90-eb5d0abd679c",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(\n",
    "    adata,\n",
    "    color=[\"leiden_res_0.02\", \"leiden_res_0.50\", \"leiden_res_2.00\"],\n",
    "    legend_loc=\"on data\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3341b3e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.umap(\n",
    "    adata,\n",
    "    color=[\"leiden_res_0.02\", \"leiden_res_0.50\", \"leiden_res_2.00\"],\n",
    "    legend_loc=\"on data\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5ca6a0e-c46c-4b7d-859f-e9c971f05d5a",
   "metadata": {},
   "source": [
    "Though UMAPs should not be over-interpreted, here we can already see that in the highest resolution our data is over-clustered, while the lowest resolution is likely grouping cells which belong to distinct cell identities."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7cac979-c45d-4639-ac0a-18f85aa383e1",
   "metadata": {},
   "source": [
    "### Marker gene set"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbc891a7-e0e8-45e8-97d9-46eb74f55cde",
   "metadata": {},
   "source": [
    "Let's define a set of marker genes for the main cell types that we expect to see in this dataset. These were adapted from [Single Cell Best Practices annotation chapter](https://www.sc-best-practices.org/cellular_structure/annotation.html), for a more detailed overview and best practices in cell type annotation, we refer the user to it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4b1f83e-c4e2-48e4-954f-72abac5b5b21",
   "metadata": {},
   "outputs": [],
   "source": [
    "marker_genes = {\n",
    "    \"CD14+ Mono\": [\"FCN1\", \"CD14\"],\n",
    "    \"CD16+ Mono\": [\"TCF7L2\", \"FCGR3A\", \"LYN\"],\n",
    "    # Note: DMXL2 should be negative\n",
    "    \"cDC2\": [\"CST3\", \"COTL1\", \"LYZ\", \"DMXL2\", \"CLEC10A\", \"FCER1A\"],\n",
    "    \"Erythroblast\": [\"MKI67\", \"HBA1\", \"HBB\"],\n",
    "    # Note HBM and GYPA are negative markers\n",
    "    \"Proerythroblast\": [\"CDK6\", \"SYNGR1\", \"HBM\", \"GYPA\"],\n",
    "    \"NK\": [\"GNLY\", \"NKG7\", \"CD247\", \"FCER1G\", \"TYROBP\", \"KLRG1\", \"FCGR3A\"],\n",
    "    \"ILC\": [\"ID2\", \"PLCG2\", \"GNLY\", \"SYNE1\"],\n",
    "    \"Naive CD20+ B\": [\"MS4A1\", \"IL4R\", \"IGHD\", \"FCRL1\", \"IGHM\"],\n",
    "    # Note IGHD and IGHM are negative markers\n",
    "    \"B cells\": [\n",
    "        \"MS4A1\",\n",
    "        \"ITGB1\",\n",
    "        \"COL4A4\",\n",
    "        \"PRDM1\",\n",
    "        \"IRF4\",\n",
    "        \"PAX5\",\n",
    "        \"BCL11A\",\n",
    "        \"BLK\",\n",
    "        \"IGHD\",\n",
    "        \"IGHM\",\n",
    "    ],\n",
    "    \"Plasma cells\": [\"MZB1\", \"HSP90B1\", \"FNDC3B\", \"PRDM1\", \"IGKC\", \"JCHAIN\"],\n",
    "    # Note PAX5 is a negative marker\n",
    "    \"Plasmablast\": [\"XBP1\", \"PRDM1\", \"PAX5\"],\n",
    "    \"CD4+ T\": [\"CD4\", \"IL7R\", \"TRBC2\"],\n",
    "    \"CD8+ T\": [\"CD8A\", \"CD8B\", \"GZMK\", \"GZMA\", \"CCL5\", \"GZMB\", \"GZMH\", \"GZMA\"],\n",
    "    \"T naive\": [\"LEF1\", \"CCR7\", \"TCF7\"],\n",
    "    \"pDC\": [\"GZMB\", \"IL3RA\", \"COBLL1\", \"TCF4\"],\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96995645-6024-4663-ad45-79e7e3aaea1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.dotplot(adata, marker_genes, groupby=\"leiden_res_0.02\", standard_scale=\"var\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e90fcd55",
   "metadata": {},
   "source": [
    "There are fairly clear patterns of expression for our markers show here, which we can use to label our coarsest clustering with broad lineages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45b3f608",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obs[\"cell_type_lvl1\"] = adata.obs[\"leiden_res_0.02\"].map(\n",
    "    {\n",
    "        \"0\": \"Lymphocytes\",\n",
    "        \"1\": \"Monocytes\",\n",
    "        \"2\": \"Erythroid\",\n",
    "        \"3\": \"B Cells\",\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a49acdff-5606-4979-b08f-e6e314b782a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.dotplot(adata, marker_genes, groupby=\"leiden_res_0.50\", standard_scale=\"var\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f6b34f4-2c5f-48ba-8b59-d0838ac1bfa3",
   "metadata": {},
   "source": [
    "This seems like a resolution that suitable to distinguish most of the different cell types in our data. As such, let's try to annotate those by manually using the dotplot above, together with the UMAP of our clusters. Ideally, one would also look specifically into each cluster, and attempt to subcluster those if required."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b904fd57-ae61-4cfa-acae-82fc36631179",
   "metadata": {},
   "source": [
    "### Differentially-expressed Genes as Markers"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0578a42c-c43b-484e-8a82-3b9ed188b94c",
   "metadata": {},
   "source": [
    "Furthermore, one can also calculate marker genes per cluster and then look up whether we can link those marker genes to any known biology, such as cell types and/or states. This is typically done using simple statistical tests, such as Wilcoxon and t-test, for each cluster vs the rest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b69a65c-199f-4883-9037-1df388b06861",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Obtain cluster-specific differentially expressed genes\n",
    "sc.tl.rank_genes_groups(adata, groupby=\"leiden_res_0.50\", method=\"wilcoxon\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "243b74f5-dd51-457d-9715-b9a0ccbfff4a",
   "metadata": {},
   "source": [
    "We can then visualize the top 5 differentially-expressed genes on a dotplot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adf45dad-5f3d-465b-acd7-aa4f8bc42f65",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.rank_genes_groups_dotplot(\n",
    "    adata, groupby=\"leiden_res_0.50\", standard_scale=\"var\", n_genes=5\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee9b9e00",
   "metadata": {},
   "source": [
    "We can then use these genes to figure out what cell types we're looking at. For example, Cluster 7 is expressing [*NKG7*](https://www.genecards.org/cgi-bin/carddisp.pl?gene=NKG7&keywords=nkg7) and [*GNLY*](https://www.genecards.org/cgi-bin/carddisp.pl?gene=GNLY&keywords=GNLY), suggesting these are [NK cells](https://en.wikipedia.org/wiki/Natural_killer_cell)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2dd5c7fc",
   "metadata": {},
   "source": [
    "To create your own plots, or use a more automated approach, the differentially expressed genes can be extracted in a convenient format with {func}`scanpy.get.rank_genes_groups_df`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fbf24a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.get.rank_genes_groups_df(adata, group=\"7\").head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef390e87",
   "metadata": {},
   "outputs": [],
   "source": [
    "dc_cluster_genes = sc.get.rank_genes_groups_df(adata, group=\"7\").head(5)[\"names\"]\n",
    "sc.pl.umap(\n",
    "    adata,\n",
    "    color=[*dc_cluster_genes, \"leiden_res_0.50\"],\n",
    "    legend_loc=\"on data\",\n",
    "    frameon=False,\n",
    "    ncols=3,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "789cdf32",
   "metadata": {},
   "source": [
    "You may have noticed that the p-values found here are extremely low. This is due to the statistical test being performed considering each cell as an independent sample. For a more conservative approach you may want to consider \"pseudo-bulking\" your data by sample (*e.g.* `sc.get.aggregate(adata, by=[\"sample\", \"cell_type\"], func=\"sum\", layer=\"counts\")`) and using a more powerful differential expression tool, like [`pydeseq2`](https://pydeseq2.readthedocs.io/)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "447eca32-2e79-4723-86f7-9ac6ab234325",
   "metadata": {},
   "source": [
    "# HoloViz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a081448-3e33-454e-98c0-7b4a3c6c6d0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# adata.write(filename='adata.h5ad')\n",
    "adata = ad.read_h5ad('adata.h5ad')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48b8f510-5c40-4458-9bb7-ff971422a0cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.var_names"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8dd79245-c308-4198-a8e2-232cf368d010",
   "metadata": {},
   "source": [
    "## Notes on anndata object:\n",
    "\n",
    "adata \n",
    "- X : main data as n_obs × n_vars array = 17041 × 23427\n",
    "- obs (observation/cell vectors as pandas Series): 'sample', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'total_counts_hb', 'log1p_total_counts_hb', 'pct_counts_hb', 'n_genes', 'doublet_score', 'predicted_doublet', 'leiden', 'leiden_res_0.02', 'leiden_res_0.50', 'leiden_res_2.00', 'cell_type_lvl1'\n",
    "- var (variable/gene vectors as pandas Series): 'mt', 'ribo', 'hb', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'highly_variable_nbatches', 'highly_variable_intersection'\n",
    "- uns (unstructured meta, config, stats as dicts): 'scrublet', 'log1p', 'hvg', 'pca', 'sample_colors', 'neighbors', 'umap', 'leiden', 'leiden_colors', 'predicted_doublet_colors', 'leiden_res_0.02', 'leiden_res_0.50', 'leiden_res_2.00', 'leiden_res_0.02_colors', 'leiden_res_0.50_colors', 'leiden_res_2.00_colors', 'rank_genes_groups', 'dendrogram_leiden_res_0.50'\n",
    "- obsm (observation/cell matrices as n_obs x whatever array): 'X_pca', 'X_umap'\n",
    "- varm (variable/gene matrices as n_var x whatever array): 'PCs'\n",
    "- layers (array transforms, same size as main data): 'counts'\n",
    "- obsp (paired observation/cell annotations as n_obs x n_obs array): 'distances', 'connectivities'\n",
    "- var_names (variable/gene vectors as pandas Index): names of genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7cca5235-0013-4183-b511-9c4443784d9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.uns['leiden_res_0.50']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bfd991f-0ca3-445c-8609-22385749eb04",
   "metadata": {},
   "source": [
    "## Convert to pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8b01ad1-2353-467f-a590-49d2d38785a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import hvplot.xarray\n",
    "import hvplot.pandas\n",
    "import holoviews as hv\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from holoviews.operation.datashader import datashade, rasterize\n",
    "hv.extension('bokeh')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1e55b33-1f1e-4527-ac1c-fb2b0787aab0",
   "metadata": {},
   "outputs": [],
   "source": [
    "umap_df = pd.DataFrame(adata.obsm['X_umap'], columns=['UMAP1', 'UMAP2'])\n",
    "pca_df = pd.DataFrame(adata.obsm['X_pca'], columns=[f'PCA{1+i}' for i in range(adata.obsm['X_pca'].shape[-1])])\n",
    "\n",
    "obs_df = adata.obs.join(umap_df.set_index(adata.obs.index))\n",
    "obs_df =  obs_df.join(pca_df.set_index(adata.obs.index))\n",
    "var_df = adata.var.copy()\n",
    "\n",
    "# Extract expression data for marker genes\n",
    "sel_genes = marker_genes['CD16+ Mono'] #[\"TCF7L2\", \"FCGR3A\", \"LYN\"],\n",
    "expression_df = pd.DataFrame(\n",
    "    adata[:, sel_genes].X.toarray(), \n",
    "    columns=sel_genes, \n",
    "    index=adata.obs_names\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b5cc1f1-ab1b-45f9-bd37-ffc7e5adf290",
   "metadata": {},
   "source": [
    "## Violins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cebe7893-72fe-4ecb-8319-b09d778e1978",
   "metadata": {},
   "outputs": [],
   "source": [
    "# later try to datshade jittered scatter as overlay using .opts(jitter=)\n",
    "cols2plot = [\"n_genes_by_counts\", \"total_counts\", \"pct_counts_mt\"]\n",
    "violins = []\n",
    "for col in cols2plot:\n",
    "    violins.append(obs_df[col].hvplot.violin(width=300, ylabel='Value', title=col))\n",
    "hv.Layout(violins).opts(shared_axes=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebeb6f2e-cecb-4806-bb60-c7e7ac29f756",
   "metadata": {},
   "source": [
    "## Scatters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecf6bf76-c6e5-438e-8da7-1d93e9f44aa0",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_df.hvplot.scatter(x=\"total_counts\", y=\"n_genes_by_counts\", color=\"pct_counts_mt\", cmap='Viridis', rasterize=True, title='pct_counts_mt', aspect=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b803ba9b-e3c9-4d45-8d87-3198fe6b7fd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# bug.. cannot get legend while using rasterize or datashade\n",
    "# bug.. cannot set color_key while using "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8812ce21-4a32-49d6-8b05-11e0898dcea4",
   "metadata": {},
   "source": [
    "## PCs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a240495e-7ac2-4ce7-a5b8-b471b95c6760",
   "metadata": {},
   "outputs": [],
   "source": [
    "pca1_2_count = obs_df.hvplot.scatter(x=\"PCA1\", y=\"PCA2\", color=\"pct_counts_mt\", cmap='Viridis', title='pct_counts_mt', rasterize=True, cnorm='eq_hist', aspect=1)\n",
    "pca1_2_samp = obs_df.hvplot.scatter(x=\"PCA1\", y=\"PCA2\", by=\"sample\", legend=True, color=['blue', 'orange'], title='sample', datashade=True, cnorm='eq_hist', aspect=1)\n",
    "# color_key={'s1d1':'blue', 's1d3':'orange'},\n",
    "pca1_2_samp + pca1_2_count"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "228c9a5b-c22c-4e4a-b0c5-8d3f410b48d6",
   "metadata": {},
   "source": [
    "## UMAP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d22490a8-838e-4a7a-9dfc-953f187a93d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "umap1_2_count = obs_df.hvplot.scatter(x=\"UMAP1\", y=\"UMAP2\", color=\"pct_counts_mt\", cmap='Viridis', title='pct_counts_mt', rasterize=True, cnorm='eq_hist', aspect=1)\n",
    "umap1_2_samp = obs_df.hvplot.scatter(x=\"UMAP1\", y=\"UMAP2\", by=\"sample\", legend=True, color=['blue', 'orange'], title='sample', datashade=True, cnorm='eq_hist', aspect=1)\n",
    "umap1_2_samp + umap1_2_count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "146ced18-ebf2-4040-b318-34a0efc87133",
   "metadata": {},
   "outputs": [],
   "source": [
    "import datashader as ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f6d4c53-35f5-49a2-9e36-d025fbaf3f06",
   "metadata": {},
   "outputs": [],
   "source": [
    "points = hv.Points(obs_df, kdims=['UMAP1', 'UMAP2'], vdims=['leiden_res_0.50'])\n",
    "umap1_2_lieden_hv = rasterize(points, aggregator=ds.by('leiden_res_0.50')).opts(\n",
    "        cmap='Category20', frame_height=300, aspect=1, tools=['hover'], title='leiden, HoloViews')\n",
    "\n",
    "umap1_2_lieden_hvplot = obs_df.hvplot.scatter(x=\"UMAP1\", y=\"UMAP2\", by=\"leiden_res_0.50\", cmap='Category20', title='leiden, hvPlot', datashade=True, cnorm='eq_hist', aspect=1)\n",
    "umap1_2_lieden_hv + umap1_2_lieden_hvplot"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "007f49e5-f3ce-4b21-92e1-59beec75a9e2",
   "metadata": {},
   "source": [
    "## dotplot with dendrogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8072394-d6b6-40f2-925d-de77cf058627",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.sparse import csr_matrix\n",
    "from scipy.cluster.hierarchy import linkage, dendrogram\n",
    "from scipy.spatial.distance import pdist\n",
    "import holoviews as hv\n",
    "from holoviews import opts\n",
    "import panel as pn\n",
    "pn.extension()\n",
    "hv.extension('bokeh')\n",
    "\n",
    "def compute_dotplot_data(expression, groupby, gene_names, markers, expression_cutoff):\n",
    "    \"\"\"\n",
    "    Compute data required for creating a dot plot, handling markers as a list or a dictionary,\n",
    "    and allowing for genes to appear in multiple groups.\n",
    "\n",
    "    Parameters:\n",
    "    - expression: n_cells x n_genes matrix (could be sparse)\n",
    "    - groupby: array of cluster assignments, length n_cells\n",
    "    - gene_names: list of gene names, length n_genes\n",
    "    - markers: either a list of gene names or a dictionary with keys as group labels and values as lists of gene names\n",
    "    - expression_cutoff: value used for binarizing expression (expression > cutoff)\n",
    "\n",
    "    Returns:\n",
    "    - df: DataFrame with columns ['gene', 'cluster', 'percentage', 'mean_expression', 'gene_group']\n",
    "    \"\"\"\n",
    "    if not isinstance(expression, csr_matrix):\n",
    "        expression = csr_matrix(expression)\n",
    "\n",
    "    gene_name_to_idx = {gene: idx for idx, gene in enumerate(gene_names)}\n",
    "\n",
    "    gene_group_list = []\n",
    "\n",
    "    if isinstance(markers, dict):\n",
    "        # Flatten markers dictionary to get list of (gene, group)\n",
    "        marker_genes = []\n",
    "        for group, genes in markers.items():\n",
    "            for gene in genes:\n",
    "                if gene in gene_name_to_idx:\n",
    "                    marker_genes.append((gene, group))\n",
    "\n",
    "        # Get the list of gene indices, allowing for duplicates\n",
    "        marker_indices = [gene_name_to_idx[gene] for gene, group in marker_genes]\n",
    "        marker_gene_names = [gene for gene, group in marker_genes]\n",
    "        gene_groups = [group for gene, group in marker_genes]\n",
    "\n",
    "    elif isinstance(markers, list):\n",
    "        marker_genes = [gene for gene in markers if gene in gene_name_to_idx]\n",
    "\n",
    "        # Remove duplicates while preserving order\n",
    "        marker_genes = list(dict.fromkeys(marker_genes))\n",
    "\n",
    "        marker_indices = [gene_name_to_idx[gene] for gene in marker_genes]\n",
    "        marker_gene_names = marker_genes\n",
    "        gene_groups = [None] * len(marker_gene_names)\n",
    "\n",
    "    else:\n",
    "        raise ValueError(\"Markers must be a list or a dictionary.\")\n",
    "\n",
    "    groupby = np.array(groupby)\n",
    "    clusters = np.unique(groupby)\n",
    "\n",
    "    clusters_series = pd.Series(clusters)\n",
    "    clusters_numeric = pd.to_numeric(clusters_series, errors='coerce')\n",
    "\n",
    "    if clusters_numeric.isnull().any():\n",
    "        convert_cluster_to_numeric = False\n",
    "    else:\n",
    "        convert_cluster_to_numeric = True\n",
    "\n",
    "    results = []\n",
    "\n",
    "    for gene_idx, gene_name, gene_group in zip(marker_indices, marker_gene_names, gene_groups):\n",
    "        gene_expression = expression[:, gene_idx]\n",
    "\n",
    "        gene_expression_binarized = gene_expression.copy()\n",
    "        gene_expression_binarized.data = (gene_expression_binarized.data > expression_cutoff).astype(int)\n",
    "\n",
    "        for cluster in clusters:\n",
    "            cluster_mask = groupby == cluster\n",
    "            cluster_cell_indices = np.where(cluster_mask)[0]\n",
    "\n",
    "            n_cells_in_cluster = len(cluster_cell_indices)\n",
    "\n",
    "            X_cluster = gene_expression[cluster_cell_indices]\n",
    "            X_cluster_binarized = gene_expression_binarized[cluster_cell_indices]\n",
    "\n",
    "            # sum of binarized expression (number of cells expressing the gene)\n",
    "            expressing_cells = X_cluster_binarized.sum()\n",
    "\n",
    "            percentage = (expressing_cells / n_cells_in_cluster) * 100\n",
    "\n",
    "            # mean expression per gene over all cells in cluster\n",
    "            total_expression = X_cluster.sum()\n",
    "            mean_expression = total_expression / n_cells_in_cluster\n",
    "\n",
    "            cluster_results = pd.DataFrame({\n",
    "                'gene': [gene_name],\n",
    "                'cluster': [cluster],\n",
    "                'percentage': [percentage],\n",
    "                'mean_expression': [mean_expression],\n",
    "                'gene_group': [gene_group]\n",
    "            })\n",
    "            results.append(cluster_results)\n",
    "\n",
    "    df = pd.concat(results, ignore_index=True)\n",
    "\n",
    "    if convert_cluster_to_numeric:\n",
    "        df['cluster'] = pd.to_numeric(df['cluster'])\n",
    "\n",
    "    return df\n",
    "\n",
    "def plot_dotplot_w_dendro(df, max_point_size=40):\n",
    "    \"\"\"\n",
    "    Create a dot plot from the DataFrame output by compute_dotplot_data,\n",
    "    handling genes that appear multiple times due to being in multiple groups,\n",
    "    and adding a dendrogram to the right of the main plot.\n",
    "\n",
    "    Parameters:\n",
    "    - df: DataFrame with columns ['gene', 'cluster', 'percentage', 'mean_expression', 'gene_group']\n",
    "    - max_point_size: Maximum size for the dots representing 100% percentage.\n",
    "\n",
    "    Returns:\n",
    "    - plot: Holoviews object representing the dot plot with a dendrogram.\n",
    "    \"\"\"\n",
    "\n",
    "    # unique identifier for each gene entry\n",
    "    df = df.copy()\n",
    "    df['gene_id'] = df.apply(\n",
    "        lambda row: f\"{row['gene']} ({row['gene_group']})\" if pd.notnull(row['gene_group']) else row['gene'],\n",
    "        axis=1\n",
    "    )\n",
    "\n",
    "    # percentage to point sizes\n",
    "    df['size'] = (df['percentage'] / df['percentage'].max()) * max_point_size\n",
    "\n",
    "    # for color mapping\n",
    "    df['mean_expression_normalized'] = df['mean_expression'] / df['mean_expression'].max()\n",
    "\n",
    "    gene_ids_order = df['gene_id'].drop_duplicates().tolist()\n",
    "    df['gene_id'] = pd.Categorical(df['gene_id'], categories=gene_ids_order, ordered=True)\n",
    "    cluster_gene_matrix = df.pivot_table(index='cluster', columns='gene_id', values='mean_expression', fill_value=0)\n",
    "\n",
    "    # hierarchical clustering\n",
    "    X = cluster_gene_matrix.values\n",
    "    cluster_dist = pdist(X, metric='euclidean')\n",
    "    cluster_linkage = linkage(cluster_dist, method='average')\n",
    "\n",
    "    dendro_data = dendrogram(cluster_linkage, labels=cluster_gene_matrix.index, no_plot=True)\n",
    "    clusters_ordered = dendro_data['ivl']\n",
    "\n",
    "    # Update df['cluster'] to be categorical with the new order\n",
    "    df['cluster'] = pd.Categorical(df['cluster'], categories=clusters_ordered, ordered=True)\n",
    "\n",
    "    # Map cluster labels to positions (for y-axis)\n",
    "    cluster_positions = {cluster: pos for pos, cluster in enumerate(clusters_ordered)}\n",
    "    df['cluster_pos'] = df['cluster'].map(cluster_positions)\n",
    "\n",
    "    cluster_values = sorted(df['cluster'].unique())\n",
    "\n",
    "    min_cluster = min(cluster_values)\n",
    "    max_cluster = max(cluster_values)\n",
    "    ylim = (max_cluster + 0.5, min_cluster - 0.5)\n",
    "\n",
    "    yticks = [(cluster, str(cluster)) for cluster in cluster_values]\n",
    "\n",
    "    points = hv.Points(\n",
    "        df,\n",
    "        kdims=['gene_id', 'cluster_pos'],\n",
    "        vdims=['mean_expression_normalized', 'size', 'percentage', 'mean_expression', 'gene_group']\n",
    "    )\n",
    "\n",
    "    points = points.opts(\n",
    "        xrotation=90,\n",
    "        color='mean_expression_normalized',\n",
    "        cmap='Reds',\n",
    "        size='size',\n",
    "        line_color='black',\n",
    "        line_alpha=0.1,\n",
    "        marker='o',\n",
    "        tools=['hover'],\n",
    "        colorbar=True,\n",
    "        width=1000,\n",
    "        frame_height=400,\n",
    "        # responsive=True,\n",
    "        xlabel='Gene',\n",
    "        ylabel='Cluster',\n",
    "        invert_yaxis=False,\n",
    "        show_legend=False,\n",
    "        yticks=yticks,\n",
    "        ylim=ylim,\n",
    "    )\n",
    "\n",
    "    if 'gene_group' in df.columns:\n",
    "        gene_groups = df[['gene_id', 'gene_group']].drop_duplicates()\n",
    "        gene_groups['group_code'] = gene_groups['gene_group'].factorize()[0]\n",
    "\n",
    "        annotations_df = gene_groups[['gene_id', 'group_code', 'gene_group']].copy()\n",
    "        annotations_df['Group'] = 'Group'\n",
    "\n",
    "        annotations_df['gene_id'] = pd.Categorical(annotations_df['gene_id'], categories=gene_ids_order, ordered=True)\n",
    "        annotations_df['Group'] = pd.Categorical(annotations_df['Group'], categories=['Group'], ordered=True)\n",
    "\n",
    "        annotations_heatmap = hv.HeatMap(\n",
    "            annotations_df,\n",
    "            kdims=['gene_id', 'Group'],\n",
    "            vdims=['group_code', 'gene_group']\n",
    "        )\n",
    "\n",
    "        annotations_heatmap = annotations_heatmap.opts(\n",
    "            colorbar=False,\n",
    "            xaxis=None,\n",
    "            yaxis=None,\n",
    "            responsive=False,\n",
    "            width=600,\n",
    "            cmap='glasbey_hv',\n",
    "            tools=['hover'],\n",
    "            toolbar=None,\n",
    "            height=50,\n",
    "            show_frame=False,\n",
    "            show_grid=False,\n",
    "        )\n",
    "\n",
    "        main_plot = hv.Layout([annotations_heatmap, hv.Empty(), points]).cols(2).opts(\n",
    "            opts.Layout(shared_axes=True)\n",
    "        )\n",
    "    else:\n",
    "        main_plot = points\n",
    "\n",
    "    # Mapping leaf positions to cluster positions\n",
    "    leaf_positions = {int(leaf_id): cluster_positions[clusters_ordered[int(leaf_id)]] for leaf_id in dendro_data['leaves']}\n",
    "\n",
    "    # What's going on? trying a swap\n",
    "    dendro_paths = []\n",
    "    icoord = np.array(dendro_data['dcoord'])\n",
    "    dcoord = np.array(dendro_data['icoord'])\n",
    "\n",
    "    for xs, ys in zip(icoord, dcoord):\n",
    "        ys_new = []\n",
    "        for y in ys:\n",
    "            if y % 10 == 5.0:\n",
    "                # Leaf node\n",
    "                leaf_id = int((y - 5.0) / 10.0)\n",
    "                ys_new.append(leaf_positions[leaf_id])\n",
    "            else:\n",
    "                # Internal node\n",
    "                ys_new.append(y / max(dcoord.flatten()) * (len(clusters_ordered) - 1))\n",
    "        dendro_paths.append(np.column_stack([xs, ys_new]))\n",
    "\n",
    "    dendrogram_plot = hv.Path(dendro_paths, ['Distance', 'Cluster'])\n",
    "\n",
    "    dendrogram_plot = dendrogram_plot.opts(\n",
    "        width=200,\n",
    "        frame_height=400,\n",
    "        xlabel='',\n",
    "        # ylabel='',\n",
    "        invert_yaxis=False,\n",
    "        xaxis=None,\n",
    "        yaxis='right',\n",
    "        # yaxis=None,\n",
    "        show_frame=False,\n",
    "        # show_grid=False,\n",
    "        fontsize={'labels': '8pt'},\n",
    "        tools=['hover'],\n",
    "        yticks=yticks,\n",
    "        ylim=ylim,\n",
    "        \n",
    "    )\n",
    "\n",
    "    return (points + dendrogram_plot).opts(shared_axes=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04e6088e-b9d4-4228-b07c-143fc9a30d3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "dp_data = compute_dotplot_data(adata.X, adata.obs['leiden_res_0.50'], adata.var_names, marker_genes, 0.1)\n",
    "plot_dotplot_w_dendro(dp_data, max_point_size=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87468954-3308-424e-bf19-4bc1ac555181",
   "metadata": {},
   "source": [
    "## dotplot with annotation bar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7487b548-34b3-45cf-ab57-3578ca1e26dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_dotplot_w_bar(df, max_point_size=40):\n",
    "    \"\"\"\n",
    "    Create a dot plot from the DataFrame output by compute_dotplot_data,\n",
    "    handling genes that appear multiple times due to being in multiple groups.\n",
    "\n",
    "    Parameters:\n",
    "    - df: DataFrame with columns ['gene', 'cluster', 'percentage', 'mean_expression', 'gene_group']\n",
    "    - max_point_size: Maximum size for the dots representing 100% percentage.\n",
    "\n",
    "    Returns:\n",
    "    - plot: Holoviews object representing the dot plot.\n",
    "    \"\"\"\n",
    "    df = df.copy()\n",
    "    df['gene_id'] = df.apply(\n",
    "        lambda row: f\"{row['gene']} ({row['gene_group']})\" if pd.notnull(row['gene_group']) else row['gene'],\n",
    "        axis=1\n",
    "    )\n",
    "\n",
    "    df['size'] = (df['percentage'] / df['percentage'].max()) * max_point_size\n",
    "    df['mean_expression_normalized'] = df['mean_expression'] / df['mean_expression'].max()\n",
    "    gene_ids_order = df['gene_id'].drop_duplicates().tolist()\n",
    "    df['gene_id'] = pd.Categorical(df['gene_id'], categories=gene_ids_order, ordered=True)\n",
    "    cluster_values = sorted(df['cluster'].unique())\n",
    "    min_cluster = min(cluster_values)\n",
    "    max_cluster = max(cluster_values)\n",
    "    ylim = (max_cluster + 0.5, min_cluster - 0.5)\n",
    "    yticks = [(cluster, str(cluster)) for cluster in cluster_values]\n",
    "    points = hv.Points(\n",
    "        df,\n",
    "        kdims=['gene_id', 'cluster'],\n",
    "        vdims=['mean_expression_normalized', 'size', 'percentage', 'mean_expression', 'gene_group']\n",
    "    )\n",
    "\n",
    "    points = points.opts(\n",
    "        xrotation=90,\n",
    "        color='mean_expression_normalized',\n",
    "        cmap='Reds',\n",
    "        size='size',\n",
    "        line_color='black',\n",
    "        line_alpha=0.1,\n",
    "        marker='o',\n",
    "        tools=['hover'],\n",
    "        colorbar=True,\n",
    "        min_height=400,\n",
    "        responsive=True,\n",
    "        xlabel='Gene',\n",
    "        ylabel='Cluster',\n",
    "        invert_yaxis=False,\n",
    "        show_legend=False,\n",
    "        yticks=yticks,\n",
    "        ylim=ylim,\n",
    "    )\n",
    "\n",
    "    if 'gene_group' in df.columns:\n",
    "        gene_groups = df[['gene_id', 'gene_group']].drop_duplicates()\n",
    "        gene_groups['group_code'] = gene_groups['gene_group'].factorize()[0]\n",
    "\n",
    "        annotations_df = gene_groups[['gene_id', 'group_code', 'gene_group']].copy()\n",
    "        annotations_df['Group'] = 'Group'\n",
    "\n",
    "        annotations_df['gene_id'] = pd.Categorical(annotations_df['gene_id'], categories=gene_ids_order, ordered=True)\n",
    "        annotations_df['Group'] = pd.Categorical(annotations_df['Group'], categories=['Group'], ordered=True)\n",
    "\n",
    "        annotations_heatmap = hv.HeatMap(\n",
    "            annotations_df,\n",
    "            kdims=['gene_id', 'Group'],\n",
    "            vdims=['group_code', 'gene_group']\n",
    "        )\n",
    "\n",
    "        annotations_heatmap = annotations_heatmap.opts(\n",
    "            colorbar=False,\n",
    "            xaxis=None,\n",
    "            yaxis=None,\n",
    "            responsive=True,\n",
    "            cmap='glasbey_hv',\n",
    "            tools=['hover'],\n",
    "            toolbar=None,\n",
    "            height=50,\n",
    "            show_frame=False,\n",
    "            show_grid=False,\n",
    "        )\n",
    "\n",
    "        layout = hv.Layout([annotations_heatmap, points]).cols(1).opts(\n",
    "            opts.Layout(shared_axes=True)\n",
    "        )\n",
    "        return layout\n",
    "\n",
    "    else:\n",
    "        return points\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10afcc96-6fd9-4f3c-bf6a-6eae486fe94a",
   "metadata": {},
   "outputs": [],
   "source": [
    "dp_data = compute_dotplot_data(adata.X, adata.obs['leiden_res_0.50'], adata.var_names, marker_genes, .1)\n",
    "plot_dotplot_w_bar(dp_data, max_point_size=12)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0922f7af-1bd6-4667-8b46-6beeb0e4a0ab",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.1"
  },
  "mystnb": {
   "execution_mode": "off"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}