Project/BraunDataset_Step2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 0,
        "height": 4,
        "hidden": false,
        "row": 0,
        "width": 4
       },
       "report_default": {}
      }
     }
    },
    "tags": []
   },
   "source": [
    "# Project: Explore Braun dataset\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Upstream Steps__\n",
    "\n",
    "* Assemble adata\n",
    "* Filtered to select cortex, cerebellum, thalamus\n",
    "* Sub-sampled to 75000 cells\n",
    "\n",
    "\n",
    "__This notebook__\n",
    "\n",
    "* QC filter on cells\n",
    "* Expression filter on genes\n",
    "* Normalization and log10 transformation by Scanpy functions\n",
    "* Feature selection (HVG) by Scanpy functions\n",
    "* Dimensionality reduction\n",
    "* Batch correction by Harmony\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dataset Information\n",
    "\n",
    "## References\n",
    "\n",
    "* [Paper Reference: Comprehensive cell atlas of the first-trimester developing human brain](https://www.science.org/doi/10.1126/science.adf1226)\n",
    "* [Data and code repository: GitHub Repo](https://github.com/linnarsson-lab/developing-human-brain)\n",
    "\n",
    "\n",
    "> Short description: 26 brain specimens spanning 5 to 14 postconceptional weeks (pcw) were dissected into 111 distinct biological samples. Each sample was subjected to single-cell RNA sequencing, resulting in 1,665,937 high-quality single-cell\n",
    "\n",
    "\n",
    "__Subsample of the dataset selecting cortex, cerebellum and thalamus__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\"><b>Cell populations identified by the authors:</b> \n",
    "\n",
    "* __Erythrocyte__\n",
    "* __Fibroblast__\n",
    "* __Glioblast__\n",
    "* __Immune__\n",
    "* __Neural crest__\n",
    "* __Neuroblast__\n",
    "* __Neuron__\n",
    "* __Neuronal IPC__\n",
    "* __Oligo__\n",
    "* __Placodes__\n",
    "* __Radial glia__\n",
    "* __Vascular__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 4,
        "height": 4,
        "hidden": false,
        "row": 0,
        "width": 4
       },
       "report_default": {}
      }
     }
    },
    "toc-hr-collapsed": true
   },
   "source": [
    "# 1. Environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scanpy as sc, anndata as ad, numpy as np, pandas as pd\n",
    "import warnings\n",
    "import yaml\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import scipy.sparse as sp\n",
    "import statsmodels.api as sm\n",
    "import scanpy as sc\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.cluster import KMeans\n",
    "\n",
    "import matplotlib.pyplot \n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "from scipy import stats\n",
    "warnings.filterwarnings('ignore')\n",
    "import plotly.express as px\n",
    "import plotly.io as pio\n",
    "import itertools\n",
    "import sys\n",
    "pio.renderers.default = \"jupyterlab\"\n",
    "\n",
    "import random\n",
    "random.seed(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "homeDir = os.getenv(\"HOME\")\n",
    "\n",
    "sc.settings.verbosity = 3             # verbosity: errors (0), warnings (1), info (2), hints (3)\n",
    "sc.logging.print_header()\n",
    "\n",
    "\n",
    "# import user defined functions from the utils folder\n",
    "import matplotlib.pyplot as plt\n",
    "sys.path.insert(1, \"../2_Day2/utils\")\n",
    "\n",
    "from CleanAdata import *\n",
    "from SankeyOBS import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.3 Input file\n",
    "\n",
    "You can set parameters useful for your analysis in the cells below and then recall them during the analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "hidden": true
       },
       "report_default": {}
      }
     }
    }
   },
   "outputs": [],
   "source": [
    "path = '../DataDir/InputData/'\n",
    "Id = 'Project_FiltNormAdata.h5ad'\n",
    "input_file = path + Id"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 0,
        "height": 4,
        "hidden": false,
        "row": 17,
        "width": 4
       },
       "report_default": {}
      }
     }
    },
    "tags": [],
    "toc-hr-collapsed": true
   },
   "source": [
    "# 1. Data Load\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata = sc.read_h5ad(input_file)\n",
    "adata.var_names_make_unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Additional QCs\n",
    "\n",
    "Here we perform few additional inspections to ensure to avoid technical issues being carried over to the next steps. \n",
    "\n",
    "First we compute again the quality metrics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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",
    "\n",
    "# hemoglobin genes\n",
    "adata.var[\"hb\"] = adata.var_names.str.contains(\"^HB[^(P)]\")\n",
    "\n",
    "sc.pp.calculate_qc_metrics(\n",
    "    adata, qc_vars=[\"mt\", \"ribo\", \"hb\"], inplace=True, log1p=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Format better Donor and Week columns\n",
    "adata.obs[\"Og_Donor\"] = adata.obs[\"Og_Donor\"].str.replace(\":\",\"\")\n",
    "adata.obs[\"Og_Week\"] = adata.obs[\"Og_Week\"].astype(\"string\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.pca(adata, color=[\"Og_Donor\", \"Og_Week\",\"Og_Chemistry\",\"Og_Subregion\"], ncols=4, wspace=.5, size = 50, vmin='p1', vmax='p99')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.4 Batches inspection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotSankey(adata, covs=[\"Og_Donor\",\"Og_Subregion\",\"Og_Chemistry\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2. Multiplets removal via scDblFinder\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import anndata2ri\n",
    "import rpy2.rinterface_lib.callbacks\n",
    "import logging\n",
    "anndata2ri.activate()\n",
    "%load_ext rpy2.ipython\n",
    "rpy2.rinterface_lib.callbacks.logger.setLevel(logging.ERROR)\n",
    "\n",
    "import os\n",
    "from rpy2 import robjects\n",
    "\n",
    "# Set R_HOME and R_LIBS_USER\n",
    "# Set R_HOME and R_LIBS_USER\n",
    "# Set R_HOME and R_LIBS_USER\n",
    "os.environ['R_HOME'] = '/opt/R/4.3.1/lib/R/bin//R'\n",
    "os.environ['R_LIBS_USER'] = '/opt/R/4.3.1/lib/R/library'\n",
    "from rpy2 import robjects\n",
    "custom_lib_paths = \"/opt/R/4.3.1/lib/R/library\"\n",
    "robjects.r(f'.libPaths(c(\"{custom_lib_paths}\", .libPaths()))')\n",
    "\n",
    "\n",
    "# Clean anndata from unnecessary fields\n",
    "# Clean anndata from unnecessary fields\n",
    "# Clean anndata from unnecessary fields\n",
    "if os.path.isfile(\"./utils/DoubletCLass.tsv\"):\n",
    "    print(\"Loading precomputed dbls class\")\n",
    "    dblfinderClass = pd.read_csv(\"./utils/DoubletCLass.tsv\", sep=\"\\t\", index_col=0)\n",
    "    adata.obs['scDblFinder.classs'] = dblfinderClass\n",
    "else:\n",
    "    dblfinderClass = pd.DataFrame()\n",
    "    for sample in adata.obs[\"sample_id\"].unique():\n",
    "        print(sample)\n",
    "        sce = adata[adata.obs[\"sample_id\"] == sample].copy()\n",
    "        sce = sce.copy()\n",
    "        sce.obs = sce.obs[[\"sample_id\"]]\n",
    "        del sce.obsp\n",
    "        sce.layers[\"counts\"] = sce.layers[\"counts\"].astype(np.int32).todense()\n",
    "        sce.X = sce.layers[\"counts\"].copy()\n",
    "        sce = sce.copy()\n",
    "\n",
    "        del sce.obsm\n",
    "        del sce.varm\n",
    "        del sce.uns\n",
    "        sce.var.index.name = None\n",
    "        sce.var = sce.var[[\"EnsembleCode\"]]\n",
    "        # Run doublets detection\n",
    "        # Run doublets detection\n",
    "        # Run doublets detection\n",
    "\n",
    "        sce = anndata2ri.py2rpy(sce)\n",
    "        print(sce)\n",
    "        scDblFinder = rpy2.robjects.packages.importr('scDblFinder')\n",
    "        S4Vectors = rpy2.robjects.packages.importr('S4Vectors')\n",
    "\n",
    "        as_data_frame = rpy2.robjects.r['as.data.frame']\n",
    "        sce = scDblFinder.scDblFinder(sce)\n",
    "        dblfinderClass = pd.concat([dblfinderClass,sce.obs[[\"scDblFinder.class\"]]])\n",
    "        #SingletsList.extend(sce.obs[\"scDblFinder.class\"].tolist())\n",
    "        del sce\n",
    "    dblfinderClass[\"scDblFinder.class\"].to_csv(\"./utils/DoubletCLass.tsv\", sep=\"\\t\")\n",
    "\n",
    "\n",
    "# Save doublets info column\n",
    "#sce.obs[\"scDblFinder.class\"].to_csv(\"...\", sep=\"\\t\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obs = pd.concat([adata.obs, dblfinderClass], axis = 1).loc[adata.obs.index]\n",
    "adata = adata[adata.obs[\"scDblFinder.class\"] == \"singlet\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3 Additional filtering ?\n",
    "\n",
    "For downstream analysis we want to be more conservative, so we inspect the relationship between the percentage of mitochondrial genes and the number of genes detected in each cell to identify outliers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-------------\n",
    "\n",
    "# 7. Additional filtering\n",
    "\n",
    "### ⚠️❓ Is it worthed to filter more this time?\n",
    "### ⚠️❓ It could be a good idea to see how clusters and dimensionality reduction is driven by any quality metric\n",
    "### ⚠️❓ Try to explore how PCs are affected by metadata (covariates)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# refine filtering\n",
    "\n",
    "sc.pl.scatter(adata, x=\"pct_counts_mt\",y=\"n_genes_by_counts\", size=40, color=\"pct_counts_ribo\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.pl.pca(adata, color=[\"Og_Donor\",\"Og_Subregion\",\"Og_Chemistry\",\"cell_class\",\"Og_Week\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that cells with a low number of genes have higher percentage of mitochondrial and ribosomal genes, so we filter out those:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata = adata[adata.obs[\"n_genes_by_counts\"] >= 2000]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4.1 PCA regressors to check variance associated with covariates\n",
    "\n",
    "A useful assessment consists in understanding how much of the variability of the PCA is explained by the covariates (**\"Donor_region\",\"Auth_Batch\",\"Auth_Assay_formatted\",\"Auth_Age\",\"cell_label\"**). We can use [Ordinary Least Squares regression](\"https://www.statsmodels.org/dev/generated/statsmodels.regression.linear_model.OLS.html\") on principal component (PC) embeddings and plot the residuals for specific covariates. This will help us understand whether for a specific principal component and a specific covariates, we observe differences across the values of the covariate. This will highlight covariates and batches that may impact the PC and answer the question __\"How much technical and biological factors guide the dimensionality reduction?\"__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from statsmodels.formula.api import ols\n",
    "\n",
    "\n",
    "def plotResiduals(adata, covToTest, npcs):\n",
    "    PCRegrDict = {}\n",
    "    sns.set_theme(style=\"ticks\")\n",
    "    varianceSummary = pd.DataFrame()\n",
    "    \n",
    "    for i in range(npcs):\n",
    "        \n",
    "        for n,c in enumerate(covToTest):\n",
    "            # format the data\n",
    "            Dummies = pd.get_dummies(adata.obs[c])\n",
    "            PCRegr = (Dummies.T * adata.obsm[\"X_pca\"][:,i].T).T.melt()\n",
    "            PCRegr = PCRegr[PCRegr[\"value\"] != 0]\n",
    "            PCRegr[\"variable\"] = PCRegr[\"variable\"].astype(\"category\")\n",
    "            PCRegr[\"cov\"] = c\n",
    "            PCRegrDict[\"PC:{}_COV:{}\".format(i+1,c)] = PCRegr.copy()\n",
    "            sns.despine(offset=30)\n",
    "            PCRegrModel = pd.get_dummies(PCRegrDict[\"PC:{}_COV:{}\".format(i+1,c)], \"variable\").copy()\n",
    "            \n",
    "            # define the regression formula\n",
    "            formula = \"\".join([\"value ~ \",\" + \".join([\"C(\"+c+\")\" for c in PCRegrModel.columns[1:-1].tolist()])])\n",
    "            \n",
    "            # fit the regression\n",
    "            fit = ols(formula, data=PCRegrModel).fit() \n",
    "            # fit.rsquared_adj\n",
    "            \n",
    "            # get the residuals\n",
    "            PCRegr[\"rsquared_adj\"] = fit.rsquared_adj\n",
    "            PCRegr[\"PC\"] = i\n",
    "            varianceSummary = pd.concat([varianceSummary,PCRegr ], ignore_index=True)\n",
    "            \n",
    "    CovOrder = {i:varianceSummary.drop_duplicates([\"PC\",\"cov\"])[varianceSummary.drop_duplicates([\"PC\",\"cov\"])[\"PC\"] == i].sort_values(\"rsquared_adj\", ascending=False)[\"cov\"].tolist() for i in varianceSummary[\"PC\"].unique().tolist()}\n",
    "    \n",
    "    for i in range(npcs):\n",
    "        fig, axes = plt.subplots(1,len(covToTest), figsize=(40,4), dpi=300, sharex=False,sharey=False)\n",
    "        plt.subplots_adjust(wspace=.5)\n",
    "        \n",
    "        for n,c in enumerate(CovOrder[i]):\n",
    "            PCRegr = varianceSummary[(varianceSummary[\"PC\"] == i) & (varianceSummary[\"cov\"] == c)]\n",
    "            sns.boxplot(data=PCRegr, x=\"variable\", y=\"value\", ax = axes[n],whis=[0, 100],\n",
    "                        palette={i:adata.uns[c+\"_colors\"][adata.obs[c].cat.categories.tolist().index(i)] for i in PCRegr[\"variable\"].unique().tolist()}\n",
    "                        #order=[i for i in PlotOrder[c] if i in PCRegr[\"variable\"].unique().tolist()],\n",
    "                        #hue_order=[i for i in PlotOrder[c] if i in PCRegr[\"variable\"].unique().tolist()])\n",
    "                       )\n",
    "            sns.stripplot(data=PCRegr, x=\"variable\", y=\"value\", size=4, color=\".3\",ax = axes[n], \n",
    "                          #order=[i for i in PlotOrder[c] if i in PCRegr[\"variable\"].unique().tolist()],\n",
    "                          #hue_order=[i for i in PlotOrder[c] if i in PCRegr[\"variable\"].unique().tolist()])\n",
    "                         )\n",
    "            axes[n].title.set_text('Covariate:{}'.format(c))\n",
    "            axes[n].set_xlabel(c, fontsize=20)\n",
    "            axes[n].set_ylabel(\"PC{} embeddings\".format(i+1), fontsize=20)\n",
    "            sns.despine(offset=30)\n",
    "            fit.rsquared_adj = PCRegr[\"rsquared_adj\"].values[0]\n",
    "            axes[n].text(0.5, 0, \"rsquared_adj:{}\".format(round(fit.rsquared_adj,2)), horizontalalignment='center', verticalalignment='center',transform=axes[n].transAxes)\n",
    "            axes[n].tick_params(axis=\"x\", rotation=45)\n",
    "            plt.xticks(rotation=45)\n",
    "            fig.autofmt_xdate(rotation=45)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npcs = 6\n",
    "covToTest = [\"Og_Donor\",\"Og_Subregion\",\"Og_Chemistry\",\"cell_class\",\"Og_Week\"]\n",
    "\n",
    "plotResiduals(adata, covToTest, npcs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. This dataset entails even more dataset, metacells could be again a good idea\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ⚠️❓ How many metacells(k) would you generate this time?\n",
    "### ⚠️❓ Feel free to play with the dataset isolating cell types and brain regions you are most interested in, or removing the ones that could drive the PCA \n",
    "### ⚠️❓ If you decide to go with metacells, which **grouping_obs** would you use?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1 Data preparation "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<div style=\"padding-top: 10px; color: red; font-size: 30px;\"> ?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2 Metacells definition and aggregation\n",
    "\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NewAdataList = []\n",
    "\n",
    "import os\n",
    "os.environ[\"OMP_NUM_THREADS\"] = \"4\"\n",
    "n_neighbors = ?\n",
    "n_pcs = ?\n",
    "n_top_genes = ?\n",
    "n_clusters = ?\n",
    "grouping_obs = ?\n",
    "\n",
    "\n",
    "#######################################\n",
    "\n",
    "if not os.path.exists(\"./utils/kmeansAggregation/\"):\n",
    "    os.makedirs(\"./utils/kmeansAggregation/\")\n",
    "\n",
    "for sample in adata.obs[grouping_obs].unique().tolist():\n",
    "    print(sample)\n",
    "    adataLocal = adata[adata.obs[grouping_obs] == sample].copy()\n",
    "    if adataLocal.shape[0] <= 400:\n",
    "        continue\n",
    "    ncells = adataLocal.shape[0]\n",
    "    filenameBackup = \"./utils/kmeansAggregation/Kmeans{}.{}.{}topgenes.{}neighbs.{}pcs.{}Ks.csv\".format(ncells,sample,n_top_genes,n_neighbors,n_pcs,n_clusters )\n",
    "    \n",
    "    # Check if the kmeans was precomputed\n",
    "    if os.path.isfile(filenameBackup):\n",
    "        print(\"Retrieving precomputed kmeans classes \\n\")\n",
    "        KmeansOBS = pd.read_csv(filenameBackup, sep=\"\\t\", index_col=0)\n",
    "        adataLocal.obs['kmeans_clusters'] = KmeansOBS\n",
    "    else:\n",
    "        #If not run k-means metacells aggregation\n",
    "        sc.pp.highly_variable_genes(adataLocal, n_top_genes=n_top_genes, flavor=\"seurat\")\n",
    "        sc.tl.pca(adataLocal)\n",
    "        # Step 1: Compute the KNN graph\n",
    "        sc.pp.neighbors(adataLocal, n_neighbors = n_neighbors , transformer='pynndescent', n_pcs=n_pcs, metric=\"euclidean\")  # Adjust n_neighbors as needed\n",
    "        # Step 2: Extract the connectivity matrix from AnnData\n",
    "        connectivities = adataLocal.obsp['distances'].toarray()\n",
    "        \n",
    "        # Step 3: Apply K-Means clustering with a fixed number of clusters\n",
    "        print(\"Computing kmeans\")\n",
    "        adataLocal.obs['kmeans_clusters'] = KMeans(n_clusters=n_clusters, random_state=0).fit_predict(connectivities)\n",
    "        print(\"Storing kmeans classes \\n\")\n",
    "        adataLocal.obs[\"kmeans_clusters\"].to_csv(filenameBackup, sep=\"\\t\")\n",
    "\n",
    "    \n",
    "    adataLocal.X = adataLocal.layers[\"counts\"].copy()\n",
    "\n",
    "    # Normalize counts\n",
    "    sc.pp.normalize_total(adataLocal, target_sum=2e4)\n",
    "\n",
    "    # Create combined AnnData objects efficiently and add to NewAdataList\n",
    "    # Step 2: Create dummy variables for clusters\n",
    "    dummy_clusters = pd.get_dummies(adataLocal.obs['kmeans_clusters'])\n",
    "\n",
    "    # Step 3: Dot product for mean aggregation of expression values\n",
    "    # Each column in `cluster_aggregated_X` represents mean expression for a cluster\n",
    "    X_dense = adataLocal.X.A if hasattr(adataLocal.X, \"A\") else adataLocal.X\n",
    "    cluster_aggregated_X = dummy_clusters.T.dot(X_dense)\n",
    "\n",
    "    cluster_aggregated_median_X = np.zeros((dummy_clusters.shape[1], X_dense.shape[1]))\n",
    "    for cluster_idx in range(dummy_clusters.shape[1]):\n",
    "        # Select cells belonging to the current cluster\n",
    "        cluster_cells = X_dense[dummy_clusters.iloc[:, cluster_idx].values == 1]\n",
    "        # Compute the median across cells within the cluster\n",
    "        cluster_aggregated_median_X[cluster_idx, :] = np.median(cluster_cells, axis=0)\n",
    "\n",
    "    # Convert to AnnData structure\n",
    "    adataAggregated = ad.AnnData(X=cluster_aggregated_X)\n",
    "    adataAggregated.layers[\"median\"] = cluster_aggregated_median_X # we save an additional layer having as aggregated value the median of the gene expression\n",
    "    adataAggregated.var_names = adataLocal.var_names\n",
    "    adataAggregated.obs['kmeans_clusters'] = dummy_clusters.columns\n",
    "\n",
    "    # Step 4: Aggregating labels and metadata\n",
    "    # Get the most common cell label for each cluster\n",
    "    adataAggregated.obs['AggregatedClass'] = (\n",
    "        adataLocal.obs.groupby('kmeans_clusters')['cell_class']\n",
    "        .agg(lambda x: x.value_counts().idxmax())\n",
    "        .reindex(dummy_clusters.columns)\n",
    "        .values\n",
    "    )\n",
    "    adataAggregated.obs['AggregatedLabel'] = (\n",
    "        adataLocal.obs.groupby('kmeans_clusters')['cell_label']\n",
    "        .agg(lambda x: x.value_counts().idxmax())\n",
    "        .reindex(dummy_clusters.columns)\n",
    "        .values\n",
    "    )\n",
    "    adataAggregated.obs_names = list(sample+\"_\"+adataAggregated.obs[\"kmeans_clusters\"].astype(str))\n",
    "    # Assign metadata fields with identical values across each cluster\n",
    "    for obsMD in [col for col in adataLocal.obs.columns if len(adataLocal.obs[col].unique()) == 1]:\n",
    "        adataAggregated.obs[obsMD] = adataLocal.obs[obsMD].iloc[0]\n",
    "\n",
    "    # Add the aggregated AnnData to NewAdataList\n",
    "    NewAdataList.append(adataAggregated)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save the new aggregated anndata\n",
    "from scipy import sparse\n",
    "CombinedAdata = ad.concat(NewAdataList)\n",
    "# Conevert to sparse matrix\n",
    "CombinedAdata.X = sparse.csr_matrix(CombinedAdata.X)\n",
    "CombinedAdata.layers[\"median\"] = sparse.csr_matrix(CombinedAdata.layers[\"median\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(adata.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(CombinedAdata.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Project continuation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"padding-top: 10px; color: green; font-weight: bold; font-size: 30px;\"> \n",
    "If you want to keep exploring the dataset by employing what we did on Day 2. Additional methods we haven't explored are welcomed!\n",
    "</div>\n"
   ]
  }
 ],
 "metadata": {
  "extensions": {
   "jupyter_dashboards": {
    "activeView": "grid_default",
    "version": 1,
    "views": {
     "grid_default": {
      "cellMargin": 10,
      "defaultCellHeight": 20,
      "maxColumns": 12,
      "name": "grid",
      "type": "grid"
     },
     "report_default": {
      "name": "report",
      "type": "report"
     }
    }
   }
  },
  "kernelspec": {
   "display_name": "Python 3",
   "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.10.12"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}