{
  "cells": [
    {
      "cell_type": "raw",
      "metadata": {},
      "source": [
        "---\n",
        "description: Grouping individual cells with similar gene expression\n",
        "  profiles to uncover distinct cell populations and their functional\n",
        "  characteristics.\n",
        "subtitle:  Scanpy Toolkit\n",
        "title:  Clustering\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div>\n",
        "\n",
        "> **Note**\n",
        ">\n",
        "> Code chunks run Python commands unless it starts with `%%bash`, in\n",
        "> which case, those chunks run shell commands.\n",
        "\n",
        "</div>\n",
        "\n",
        "In this tutorial we will continue the analysis of the integrated\n",
        "dataset. We will use the scanpy enbedding to perform the clustering\n",
        "using graph community detection algorithms.\n",
        "\n",
        "Let's first load all necessary libraries and also the integrated dataset\n",
        "from the previous step."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "import numpy as np\n",
        "import pandas as pd\n",
        "import scanpy as sc\n",
        "import matplotlib.pyplot as plt\n",
        "import warnings\n",
        "import os\n",
        "import urllib.request\n",
        "\n",
        "warnings.simplefilter(action=\"ignore\", category=Warning)\n",
        "\n",
        "# verbosity: errors (0), warnings (1), info (2), hints (3)\n",
        "sc.settings.verbosity = 3\n",
        "sc.settings.set_figure_params(dpi=80)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# download pre-computed data if missing or long compute\n",
        "fetch_data = True\n",
        "\n",
        "# url for source and intermediate data\n",
        "path_data = \"https://export.uppmax.uu.se/naiss2023-23-3/workshops/workshop-scrnaseq\"\n",
        "\n",
        "path_results = \"data/covid/results\"\n",
        "if not os.path.exists(path_results):\n",
        "    os.makedirs(path_results, exist_ok=True)\n",
        "\n",
        "path_file = \"data/covid/results/scanpy_covid_qc_dr_scanorama.h5ad\"\n",
        "if fetch_data and not os.path.exists(path_file):\n",
        "    urllib.request.urlretrieve(os.path.join(\n",
        "        path_data, 'covid/results/scanpy_covid_qc_dr_scanorama.h5ad'), path_file)\n",
        "\n",
        "adata = sc.read_h5ad(path_file)\n",
        "adata"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Graph clustering\n",
        "\n",
        "The procedure of clustering on a Graph can be generalized as 3 main\n",
        "steps:\\\n",
        "- Build a kNN graph from the data.\\\n",
        "- Prune spurious connections from kNN graph (optional step). This is a\n",
        "SNN graph.\\\n",
        "- Find groups of cells that maximizes the connections within the group\n",
        "compared other groups.\n",
        "\n",
        "If you recall from the integration, we already constructed a knn graph\n",
        "before running UMAP. Hence we do not need to do it again, and can run\n",
        "the community detection right away.\n",
        "\n",
        "The modularity optimization algoritm in Scanpy are *Leiden* and\n",
        "*Louvain*. Lets test both and see how they compare.\n",
        "\n",
        "### Leiden"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "sc.tl.leiden(adata, key_added = \"leiden_1.0\") # default resolution in 1.0\n",
        "sc.tl.leiden(adata, resolution = 0.6, key_added = \"leiden_0.6\")\n",
        "sc.tl.leiden(adata, resolution = 0.4, key_added = \"leiden_0.4\")\n",
        "sc.tl.leiden(adata, resolution = 1.4, key_added = \"leiden_1.4\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the clusters, as you can see, with increased resolution, we get\n",
        "higher granularity in the clustering."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "sc.pl.umap(adata, color=['leiden_0.4', 'leiden_0.6', 'leiden_1.0','leiden_1.4'])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Once we have done clustering, the relationships between clusters can be\n",
        "calculated as correlation in PCA space and we also visualize some of the\n",
        "marker genes that we used in the Dim Reduction lab onto the clusters."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "sc.tl.dendrogram(adata, groupby = \"leiden_0.6\")\n",
        "sc.pl.dendrogram(adata, groupby = \"leiden_0.6\")\n",
        "\n",
        "genes  = [\"CD3E\", \"CD4\", \"CD8A\", \"GNLY\",\"NKG7\", \"MS4A1\",\"FCGR3A\",\"CD14\",\"LYZ\",\"CST3\",\"MS4A7\",\"FCGR1A\"]\n",
        "sc.pl.dotplot(adata, genes, groupby='leiden_0.6', dendrogram=True)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Louvain"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "sc.tl.louvain(adata, key_added = \"louvain_1.0\") # default resolution in 1.0\n",
        "sc.tl.louvain(adata, resolution = 0.6, key_added = \"louvain_0.6\")\n",
        "sc.tl.louvain(adata, resolution = 0.4, key_added = \"louvain_0.4\")\n",
        "sc.tl.louvain(adata, resolution = 1.4, key_added = \"louvain_1.4\")\n",
        "\n",
        "sc.pl.umap(adata, color=['louvain_0.4', 'louvain_0.6', 'louvain_1.0','louvain_1.4'])\n",
        "\n",
        "sc.tl.dendrogram(adata, groupby = \"louvain_0.6\")\n",
        "sc.pl.dendrogram(adata, groupby = \"louvain_0.6\")\n",
        "\n",
        "genes  = [\"CD3E\", \"CD4\", \"CD8A\", \"GNLY\",\"NKG7\", \"MS4A1\",\"FCGR3A\",\"CD14\",\"LYZ\",\"CST3\",\"MS4A7\",\"FCGR1A\"]\n",
        "\n",
        "sc.pl.dotplot(adata, genes, groupby='louvain_0.6', dendrogram=True)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## K-means clustering\n",
        "\n",
        "K-means is a generic clustering algorithm that has been used in many\n",
        "application areas. In R, it can be applied via the `kmeans()` function.\n",
        "Typically, it is applied to a reduced dimension representation of the\n",
        "expression data (most often PCA, because of the interpretability of the\n",
        "low-dimensional distances). We need to define the number of clusters in\n",
        "advance. Since the results depend on the initialization of the cluster\n",
        "centers, it is typically recommended to run K-means with multiple\n",
        "starting configurations (via the `nstart` argument)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "from sklearn.cluster import KMeans\n",
        "from sklearn.metrics import adjusted_rand_score\n",
        "\n",
        "# extract pca coordinates\n",
        "X_pca = adata.obsm['Scanorama'] \n",
        "\n",
        "# kmeans with k=5\n",
        "kmeans = KMeans(n_clusters=5, random_state=0).fit(X_pca) \n",
        "adata.obs['kmeans5'] = kmeans.labels_.astype(str)\n",
        "\n",
        "# kmeans with k=10\n",
        "kmeans = KMeans(n_clusters=10, random_state=0).fit(X_pca) \n",
        "adata.obs['kmeans10'] = kmeans.labels_.astype(str)\n",
        "\n",
        "# kmeans with k=15\n",
        "kmeans = KMeans(n_clusters=15, random_state=0).fit(X_pca)\n",
        "adata.obs['kmeans15'] = kmeans.labels_.astype(str)\n",
        "\n",
        "sc.pl.umap(adata, color=['kmeans5', 'kmeans10', 'kmeans15'])\n",
        "\n",
        "adata.obsm"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Hierarchical clustering\n",
        "\n",
        "Hierarchical clustering is another generic form of clustering that can\n",
        "be applied also to scRNA-seq data. As K-means, it is typically applied\n",
        "to a reduced dimension representation of the data. Hierarchical\n",
        "clustering returns an entire hierarchy of partitionings (a dendrogram)\n",
        "that can be cut at different levels. Hierarchical clustering is done in\n",
        "these steps:\n",
        "\n",
        "1.  Define the distances between samples. The most common are Euclidean\n",
        "    distance (a.k.a. straight line between two points) or correlation\n",
        "    coefficients.\n",
        "2.  Define a measure of distances between clusters, called *linkage*\n",
        "    criteria. It can for example be average distances between clusters.\n",
        "    Commonly used methods are `single`, `complete`, `average`, `median`,\n",
        "    `centroid` and `ward`.\n",
        "3.  Define the dendrogram among all samples using **Bottom-up** or\n",
        "    **Top-down** approach. **Bottom-up** is where samples start with\n",
        "    their own cluster which end up merged pair-by-pair until only one\n",
        "    cluster is left. **Top-down** is where samples start all in the same\n",
        "    cluster that end up being split by 2 until each sample has its own\n",
        "    cluster.\n",
        "\n",
        "As you might have realized, correlation is not a method implemented in\n",
        "the `dist()` function. However, we can create our own distances and\n",
        "transform them to a distance object. We can first compute sample\n",
        "correlations using the `cor` function.\\\n",
        "As you already know, correlation range from -1 to 1, where 1 indicates\n",
        "that two samples are closest, -1 indicates that two samples are the\n",
        "furthest and 0 is somewhat in between. This, however, creates a problem\n",
        "in defining distances because a distance of 0 indicates that two samples\n",
        "are closest, 1 indicates that two samples are the furthest and distance\n",
        "of -1 is not meaningful. We thus need to transform the correlations to a\n",
        "positive scale (a.k.a. **adjacency**):\\\n",
        "$$adj = \\frac{1- cor}{2}$$\\\n",
        "Once we transformed the correlations to a 0-1 scale, we can simply\n",
        "convert it to a distance object using `as.dist()` function. The\n",
        "transformation does not need to have a maximum of 1, but it is more\n",
        "intuitive to have it at 1, rather than at any other number.\n",
        "\n",
        "The function `AgglomerativeClustering` has the option of running with\n",
        "disntance metrics \"euclidean\", \"l1\", \"l2\", \"manhattan\", \"cosine\", or\n",
        "\"precomputed\". However, with ward linkage only euklidean distances\n",
        "works. Here we will try out euclidean distance and ward linkage\n",
        "calculated in PCA space."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "from sklearn.cluster import AgglomerativeClustering\n",
        "\n",
        "cluster = AgglomerativeClustering(n_clusters=5, linkage='ward')\n",
        "adata.obs['hclust_5'] = cluster.fit_predict(X_pca).astype(str)\n",
        "\n",
        "cluster = AgglomerativeClustering(n_clusters=10, linkage='ward')\n",
        "adata.obs['hclust_10'] = cluster.fit_predict(X_pca).astype(str)\n",
        "\n",
        "cluster = AgglomerativeClustering(n_clusters=15, linkage='ward')\n",
        "adata.obs['hclust_15'] = cluster.fit_predict(X_pca).astype(str)\n",
        "\n",
        "sc.pl.umap(adata, color=['hclust_5', 'hclust_10', 'hclust_15'])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, lets save the clustered data for further analysis."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "adata.write_h5ad('./data/covid/results/scanpy_covid_qc_dr_scanorama_cl.h5ad')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Distribution of clusters\n",
        "\n",
        "Now, we can select one of our clustering methods and compare the\n",
        "proportion of samples across the clusters.\n",
        "\n",
        "Select the \"leiden_0.6\" and plot proportion of samples per cluster and\n",
        "also proportion covid vs ctrl.\n",
        "\n",
        "Plot proportion of cells from each condition per cluster."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "tmp = pd.crosstab(adata.obs['leiden_0.6'],adata.obs['type'], normalize='index')\n",
        "tmp.plot.bar(stacked=True).legend(bbox_to_anchor=(1.4, 1), loc='upper right')\n",
        "\n",
        "tmp = pd.crosstab(adata.obs['leiden_0.6'],adata.obs['sample'], normalize='index')\n",
        "tmp.plot.bar(stacked=True).legend(bbox_to_anchor=(1.4, 1),loc='upper right')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In this case we have quite good representation of each sample in each\n",
        "cluster. But there are clearly some biases with more cells from one\n",
        "sample in some clusters and also more covid cells in some of the\n",
        "clusters.\n",
        "\n",
        "We can also plot it in the other direction, the proportion of each\n",
        "cluster per sample."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "tmp = pd.crosstab(adata.obs['sample'],adata.obs['leiden_0.6'], normalize='index')\n",
        "tmp.plot.bar(stacked=True).legend(bbox_to_anchor=(1.4, 1), loc='upper right')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div>\n",
        "\n",
        "> **Discuss**\n",
        ">\n",
        "> By now you should know how to plot different features onto your data.\n",
        "> Take the QC metrics that were calculated in the first exercise, that\n",
        "> should be stored in your data object, and plot it as violin plots per\n",
        "> cluster using the clustering method of your choice. For example, plot\n",
        "> number of UMIS, detected genes, percent mitochondrial reads. Then,\n",
        "> check carefully if there is any bias in how your data is separated by\n",
        "> quality metrics. Could it be explained biologically, or could there be\n",
        "> a technical bias there?\n",
        "\n",
        "</div>\n",
        "\n",
        "## Session info\n",
        "\n",
        "```{=html}\n",
        "<details>\n",
        "```\n",
        "```{=html}\n",
        "<summary>\n",
        "```\n",
        "Click here\n",
        "```{=html}\n",
        "</summary>\n",
        "```"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "sc.logging.print_versions()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "```{=html}\n",
        "</details>\n",
        "```"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}