Coarse_cell_typing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "824ad282-0f19-40ec-b32a-ca07dbbbdbf8",
   "metadata": {},
   "source": [
    "# Coarse Cell Typing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(warn = -1, verbose=FALSE)\n",
    "#!/usr/bin/env Rscript \n",
    "library(dplyr)\n",
    "library(Seurat)\n",
    "library(httr)\n",
    "library(readr)\n",
    "library(pheatmap)\n",
    "library(RColorBrewer)\n",
    "library(ggplot2)\n",
    "library(cowplot)\n",
    "library(patchwork)\n",
    "library(unixtools)\n",
    "library(ggrepel)\n",
    "library(repr)\n",
    "library(ggmin)\n",
    "library(harmony)\n",
    "library(SeuratWrappers)\n",
    "library(Nebulosa)\n",
    "library(ggthemes)\n",
    "library(purrr)\n",
    "library(radiant.data)\n",
    "library(presto)\n",
    "library(pryr)\n",
    "set_config(config(ssl_verifypeer = 0L))\n",
    "ulimit::memory_limit(200000)\n",
    "set.tempdir(\"/datastore/lucy/tmp/\")\n",
    "setwd(\"/datastore/lucy/CosMx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba6e9d48-6110-4c23-9115-3697983f3d8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "source(\"./R/utils.R\")\n",
    "start_upR(clusterfiles = TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "QC_harmony_pipeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RunModularityClustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RunModularityClusteringCpp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_dim_red"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0a60261-3211-4718-b804-b7294cb6927a",
   "metadata": {},
   "source": [
    "# Read-in data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a0421b6-3fcd-44a2-a22c-11dc4c574c9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "obj<-list()\n",
    "obj$counts<-readRDS(\"../data_Glasgow/Glasgow.RDS\")\n",
    "obj$metadata<-fread(\"../data_Glasgow/Glasgow.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acf214c4-2d60-48e8-a13e-b5296590675e",
   "metadata": {},
   "outputs": [],
   "source": [
    "cellgeoms<-readRDS(\"../data_Glasgow/Glasgow_cellgeoms_QC.RDS\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00058b0b-5c0b-4f15-9c60-68ac6678f228",
   "metadata": {},
   "source": [
    "## QC based on cell size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d7939d2-bcbe-420f-9b35-d0d1cf554ec1",
   "metadata": {},
   "outputs": [],
   "source": [
    "obj$metadata<-obj$metadata %>% subset(cellID %in% cellgeoms$cellID)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9d5d640-11af-4c06-952a-1ee20ef170d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "dimlist(obj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb365493-70f0-46b8-b96f-a9395becdb46",
   "metadata": {},
   "outputs": [],
   "source": [
    "head(obj$metadata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9438d23e-c3b0-48f3-b7a3-acf1571e4ac2",
   "metadata": {},
   "outputs": [],
   "source": [
    "unique(obj$metadata$Disease)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72327466-ff36-4896-90b3-7ea771f36b18",
   "metadata": {},
   "outputs": [],
   "source": [
    "all(obj$metadata$cellID == colnames(obj$counts))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "877c5c83-729b-474e-96af-bd16acfb3b1b",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Genecellmat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5c3ca1a-5052-4be5-a009-efdfdf204850",
   "metadata": {},
   "outputs": [],
   "source": [
    "dimlist(obj)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fa47a4f-9f31-40e9-b7ba-c2cb49879e3d",
   "metadata": {},
   "source": [
    "## QC - I loaded pre-QCed data: \n",
    "+ ngenes>=20 & ncounts>=30"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08dd6905-8caa-449e-a8f3-76c42e4e0626",
   "metadata": {},
   "source": [
    "## delete negprobs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c868717-b79a-4f45-9177-c7f2af715dbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "posPrbs<-rownames(obj$counts)[!grepl(\"NegPrb\", rownames(obj$counts))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0737616f-a7bd-4ddd-8350-0a8cd3e4134c",
   "metadata": {},
   "outputs": [],
   "source": [
    "length(posPrbs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "076c9b5b-703e-4551-8b9e-82519492ff54",
   "metadata": {},
   "outputs": [],
   "source": [
    "dim(obj$counts[posPrbs, ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6207216-49b0-4d6a-ba4f-2a1a213d4956",
   "metadata": {},
   "outputs": [],
   "source": [
    "obj$counts<-obj$counts[posPrbs, ]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80687326-6cb1-4b2b-8acb-f3091426b1f3",
   "metadata": {},
   "source": [
    "## ngenes per cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b549401-801b-4c86-8215-aa50bfb35708",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(5, 10)\n",
    "plot_genes_counts(obj$counts, gene_int = 20, count_int = 30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16b9651e-99c2-4a58-acaa-d34869c68eba",
   "metadata": {},
   "outputs": [],
   "source": [
    "range(colSums(obj$counts > 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e6a7c84-b0fa-41ae-994a-b1074dedc8a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "range(colSums(obj$counts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bf01e2f-3060-4751-b286-b04109f0d9a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "obj$counts[1:5, 1:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d537d3c2-c20a-4804-af73-538af28ffaa4",
   "metadata": {},
   "source": [
    "# Harmonize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff7ee2c0-84de-46c1-bb93-21b856abaf8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "norm_value<-median(colSums(obj$counts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b43d538-fe14-483b-9117-cda63105918c",
   "metadata": {},
   "outputs": [],
   "source": [
    "norm_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a43f621c-513a-486e-860c-db6e6aa66faa",
   "metadata": {},
   "outputs": [],
   "source": [
    "obj$metadata<-obj$metadata %>% \n",
    "    unite(\"SampleFOV\", c(\"SampleID\", \"FOV\"), remove = FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "060d1719-5b96-4633-ae51-b562c5ecc941",
   "metadata": {},
   "outputs": [],
   "source": [
    "head(obj$metadata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cdd218d-aaaf-47cc-80f5-f7f8583ff821",
   "metadata": {},
   "outputs": [],
   "source": [
    "batch<-c('SampleID', 'SampleFOV')\n",
    "cluster_res<-c(1.5, 2.5, 3.1)\n",
    "theta_harmony<-c(0, 0) \n",
    "sigma_harmony<-0.25\n",
    "system.time({\n",
    "    objH<-QC_harmony_pipeline(\n",
    "        obj, \n",
    "        ngenes_threshold = 20, \n",
    "        ncounts_threshold = 30, \n",
    "        norm_value, \n",
    "        do_cluster_after = TRUE,\n",
    "        resolution_clustering = cluster_res, \n",
    "        clustering_ncores = 8,\n",
    "        vars_use = batch,\n",
    "        theta = theta_harmony,\n",
    "        sigma = sigma_harmony,\n",
    "        max.iter.harmony = 15,\n",
    "        max.iter.cluster = 100,\n",
    "        return_object = TRUE\n",
    "    )\n",
    "\n",
    "    objH$sigma_harmony<-sigma_harmony\n",
    "    objH$vars_use<-batch\n",
    "    objH$theta_harmony<-theta_harmony    \n",
    "\n",
    "    })\n",
    "\n",
    "saveRDS(objH, \"../data_Glasgow/cache/CosMxcoarseGrainharmonyObj.RDS\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a80da6fe-bde1-42b7-a259-e282638a40e2",
   "metadata": {},
   "source": [
    "## load"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d2fa4b2-5dc3-4780-aa04-4bd0ca882ba3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# objH<-readRDS( \"../data_Glasgow/cache/CosMxcoarseGrainharmonyObj.RDS\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b803fbf-8893-4ef2-945e-4d1940352e99",
   "metadata": {},
   "outputs": [],
   "source": [
    "obj$metadata<-obj$Humap$embedding %>% \n",
    "    as.data.frame %>% \n",
    "    cbind(obj$Humap$clusters) %>% \n",
    "    rownames_to_column('cellID') %>% \n",
    "    left_join(obj$metadata, by = 'cellID')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87c593e5-bd63-469a-a359-44d202f578d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "head(obj$metadata)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65da6ff2-ed47-49bb-91b4-028521b4034d",
   "metadata": {},
   "source": [
    "## Figures"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c70423b1-6ee9-4fb3-8597-d9a892bbc387",
   "metadata": {},
   "source": [
    "### PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3665acc7-b8c6-4981-98d3-24f8aa652cc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "pca_before<-plot_dim_red(\n",
    "    dim_red_embeddings = objH$pca$embeddings,\n",
    "    clusters = NULL,\n",
    "    metadata = objH$metadata, \n",
    "    cell_id_colname = \"cellID\",\n",
    "    color_by = \"SampleID\",\n",
    "    plot_title = \"PCA before Harmony\",\n",
    "    dim_red_type = \"PCA\",\n",
    "    size_points = 0.2\n",
    ")\n",
    "pca_after<-plot_dim_red(\n",
    "    dim_red_embeddings = objH$H,\n",
    "    clusters = NULL,\n",
    "    metadata = objH$metadata, \n",
    "    cell_id_colname = \"cellID\",\n",
    "    color_by = \"SampleID\",\n",
    "    plot_title = \"PCA after Harmony\",\n",
    "    dim_red_type = \"PCA\",\n",
    "    size_points = 0.2\n",
    ")\n",
    "fig.size(5, 10)\n",
    "pca_before + theme(legend.position = \"none\") | \n",
    "pca_after + theme(legend.position = \"none\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e2041af-5a02-4813-96a0-9f2c78e9a63b",
   "metadata": {},
   "source": [
    "### UMAP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d301524-69d0-4d61-86ae-cea6ad6fc23b",
   "metadata": {},
   "outputs": [],
   "source": [
    "umap_before<-plot_dim_red(\n",
    "    dim_red_embeddings = objH$umap$embedding,\n",
    "    clusters = NULL,\n",
    "    metadata = objH$metadata, \n",
    "    cell_id_colname = \"cellID\",\n",
    "    color_by = \"SampleID\",\n",
    "    plot_title = \"UMAP before Harmony\",\n",
    "    dim_red_type = \"UMAP\",\n",
    "    size_points = 0.2\n",
    ")\n",
    "umap_after<-plot_dim_red(\n",
    "    dim_red_embeddings = objH$Humap$embedding,\n",
    "    clusters = NULL,\n",
    "    metadata = objH$metadata, \n",
    "    cell_id_colname = \"cellID\",\n",
    "    color_by = \"SampleID\",\n",
    "    plot_title = \"UMAP after Harmony\",\n",
    "    dim_red_type = \"UMAP\",\n",
    "    size_points = 0.2\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d36351f8-0f89-404a-bf5b-a9efe703ae1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(5, 15)\n",
    "umap_before | \n",
    "umap_after "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5279c3fc-ee2f-43c0-a340-a7b154afb7ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "umap_before<-plot_dim_red(\n",
    "    dim_red_embeddings = objH$umap$embedding,\n",
    "    clusters = NULL,\n",
    "    metadata = objH$metadata, \n",
    "    cell_id_colname = \"cellID\",\n",
    "    color_by = \"SampleFOV\",\n",
    "    plot_title = \"UMAP before Harmony\",\n",
    "    dim_red_type = \"UMAP\",\n",
    "    size_points = 0.2\n",
    ")\n",
    "umap_after<-plot_dim_red(\n",
    "    dim_red_embeddings = objH$Humap$embedding,\n",
    "    clusters = NULL,\n",
    "    metadata = objH$metadata, \n",
    "    cell_id_colname = \"cellID\",\n",
    "    color_by = \"SampleFOV\",\n",
    "    plot_title = \"UMAP after Harmony\",\n",
    "    dim_red_type = \"UMAP\",\n",
    "    size_points = 0.2\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1915962f-679f-4372-b304-7895a39e1206",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(5, 15)\n",
    "\n",
    "umap_after "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ebbe3d1-12d0-4f9e-8687-2b88fcf6a389",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(5, 8)\n",
    "umap_after<-plot_dim_red(\n",
    "    dim_red_embeddings = objH$Humap$embedding,\n",
    "    clusters = objH$Humap$clusters,\n",
    "    metadata = objH$metadata, \n",
    "    cell_id_colname = \"cellID\",\n",
    "    color_by = \"Clust1.5\",\n",
    "    plot_title = \"UMAP after Harmony\",\n",
    "    dim_red_type = \"UMAP\",\n",
    "    size_points = 0.2\n",
    ")\n",
    "umap_after"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c545d496-d415-422b-8bdb-312e07a5c908",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(7, 12)\n",
    "umap_after<-plot_dim_red(\n",
    "    dim_red_embeddings = objH$Humap$embedding,\n",
    "    clusters = objH$Humap$clusters,\n",
    "    metadata = objH$metadata, \n",
    "    cell_id_colname = \"cellID\",\n",
    "    color_by = \"Clust3.1\",\n",
    "    plot_title = \"UMAP after Harmony\",\n",
    "    dim_red_type = \"UMAP\",\n",
    "    size_points = 0.2\n",
    ")\n",
    "umap_after"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1897d7ef-2b5f-4f14-a267-63fe6e129a96",
   "metadata": {},
   "source": [
    "# Find markers"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e568071-fbbe-4229-a1fe-99a104e9e873",
   "metadata": {},
   "source": [
    "## load"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "08cd707e-336f-44d2-b59d-5443b2aa288f",
   "metadata": {},
   "source": [
    "## Clust3.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca50c77c-f353-4645-a2a4-61a786020ef0",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(5, 12)\n",
    "umap_after<-plot_dim_red(\n",
    "    dim_red_embeddings = objH$Humap$embedding,\n",
    "    clusters = objH$Humap$clusters,\n",
    "    metadata = objH$metadata, \n",
    "    cell_id_colname = \"cellID\",\n",
    "    color_by = \"Clust3.1\",\n",
    "    plot_title = \"UMAP after Harmony\",\n",
    "    dim_red_type = \"UMAP\",\n",
    "    size_points = 0.2,\n",
    "    plot_labels = FALSE\n",
    ")\n",
    "umap_after"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "172c750c-5065-409b-b405-dc5c513598ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(10, 20)\n",
    "umap_after + facet_wrap(~ color_col)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a52c056-5970-455e-aa66-6fcfd22a3122",
   "metadata": {},
   "outputs": [],
   "source": [
    "formula_glmer_input<-createFormula(\"y\", \"1 + (1|Clust3.1)\", c(\"(1|SampleFOV/Clust3.1)\", \"(1|SampleID/Clust3.1)\", \"offset(logUMI)\"))\n",
    "\n",
    "objH$markers_Clust3.1<-find_glmer_markers(\n",
    "    objH$counts,\n",
    "    objH$metadata %>% cbind(objH$Humap$clusters),\n",
    "    formula_glmer_input, \n",
    "    \"Clust3.1\", \n",
    "    c(\"SampleID\", \"FOV\", \"Clust3.1\") \n",
    ")\n",
    "saveRDS(objH, \"../data_Glasgow/cache/CosMxcoarseGrainharmonyObj_markers.RDS\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdd4e45a-9ea1-42a5-965f-7f704ed07d98",
   "metadata": {},
   "source": [
    "### Viz clusters at 3.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57803dfd-8999-4636-b01c-bc1e037e4a8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.matrix.max.cols=50, repr.matrix.max.rows=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc9b47c3-773b-4b7f-bb5a-e333756675b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "objH$markers_Clust3.1 %>% \n",
    "    mutate(contrast = as.numeric(contrast)) %>% \n",
    "    group_by(contrast) %>% \n",
    "    filter(pvalue < 1e-02) %>% \n",
    "    filter(logFC > 0.5) %>% \n",
    "    top_n(n = 20, wt = logFC) %>% \n",
    "    mutate(rank = rank(-logFC)) %>% \n",
    "    ungroup %>% \n",
    "    # dplyr::mutate(feature_logFC = paste0(feature, \"_\", logFC)) %>% \n",
    "    dplyr::select(contrast, feature, rank) %>% \n",
    "    spread(contrast, feature, fill = NA)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d6a5408-909a-46ce-a6e2-bfd7e8804b86",
   "metadata": {},
   "source": [
    "## Clust1.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a99d6d2d-a729-4076-9f24-bfbfcbe05e33",
   "metadata": {},
   "outputs": [],
   "source": [
    "formula_glmer_input<-createFormula(\"y\", \"1 + (1|Clust1.5)\", c(\"(1|SampleFOV/Clust1.5)\", \"(1|SampleID/Clust1.5)\", \"offset(logUMI)\"))\n",
    "\n",
    "objH$markers_Clust1.5<-find_glmer_markers(\n",
    "    objH$counts,\n",
    "    objH$metadata %>% cbind(objH$Humap$clusters),\n",
    "    formula_glmer_input, \n",
    "    \"Clust1.5\", \n",
    "    c(\"SampleID\", \"FOV\", \"Clust1.5\") \n",
    ")\n",
    "saveRDS(objH, \"../data_Glasgow/cache/CosMxcoarseGrainharmonyObj_markers.RDS\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objH <- readRDS(\"./cache/CosMxcoarseGrainharmonyObj_markers.RDS\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60f23b32-291c-4c86-b8e5-bbf523801a60",
   "metadata": {},
   "source": [
    "### Viz clusters at 1.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2068c896-9d58-4616-be13-20e567224b65",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(5, 8)\n",
    "umap_after<-plot_dim_red(\n",
    "    dim_red_embeddings = objH$Humap$embedding,\n",
    "    clusters = objH$Humap$clusters,\n",
    "    metadata = objH$metadata, \n",
    "    cell_id_colname = \"cellID\",\n",
    "    color_by = \"Clust1.5\",\n",
    "    plot_title = \"UMAP after Harmony\",\n",
    "    dim_red_type = \"UMAP\",\n",
    "    size_points = 0.2,\n",
    "    plot_labels = FALSE\n",
    ")\n",
    "umap_after"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objH$metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d08cc117-5e6e-4ba8-ae23-ebe8bdf699bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(10, 15)\n",
    "umap_after + facet_wrap(~ color_col)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5528e5e-f788-4278-a02c-13b0eee35989",
   "metadata": {},
   "outputs": [],
   "source": [
    "objH$markers_Clust1.5 %>% \n",
    "    mutate(contrast = as.numeric(contrast)) %>% \n",
    "    group_by(contrast) %>% \n",
    "    filter(pvalue < 1e-02) %>% \n",
    "    filter(logFC > 0.5) %>% \n",
    "    top_n(n = 20, wt = logFC) %>% \n",
    "    mutate(rank = rank(-logFC)) %>% \n",
    "    ungroup %>% \n",
    "    # dplyr::mutate(feature_logFC = paste0(feature, \"_\", logFC)) %>% \n",
    "    dplyr::select(contrast, feature, rank) %>% \n",
    "    spread(contrast, feature, fill = NA)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38afe8bb-0b4d-418e-8a17-d61af117c781",
   "metadata": {},
   "source": [
    "### annotate at Clust1.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba649ddf-8d31-4e4a-8efb-0806b479e2d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of cells in each cluster\n",
    "objH$metadata %>% \n",
    "    cbind(objH$Humap$clusters) %>% \n",
    "    with(table(Clust1.5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08186cd0-cbbd-435a-8a6d-9b79018531e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "objH$metadata<-objH$metadata %>% \n",
    "    cbind(objH$Humap$clusters) %>% \n",
    "    mutate(celltype.coarse = case_when(\n",
    "        Clust1.5 %in%(c(0)) ~ \"T cells\",\n",
    "        Clust1.5 %in%(c(1, 5, 6, 12)) ~ \"Myeloid\",\n",
    "        Clust1.5 %in%(c(7)) ~ \"DC\",\n",
    "        Clust1.5 %in%(c(2, 8, 9, 11)) ~ \"Fibroblasts\",\n",
    "        Clust1.5 %in%(c(3)) ~ \"Endothelial\",\n",
    "        Clust1.5 %in%(c(4)) ~ \"Plasma\",\n",
    "        Clust1.5 %in%(c(10)) ~ \"Mural\",\n",
    "        Clust1.5 %in%(c(13)) ~ \"B cells\",\n",
    "        Clust1.5 %in%(c(14)) ~ \"Mast\",\n",
    "        Clust1.5 %in%(c(15)) ~ \"Adipocytes\",\n",
    "        Clust1.5 %in%(c(16)) ~ \"ProliferationT\",\n",
    "        .default = Clust1.5\n",
    "        )\n",
    "    ) %>% \n",
    "    dplyr::select(- starts_with(\"Clust\"))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cc2b1fc-034c-414f-986e-629f41d90cdd",
   "metadata": {},
   "outputs": [],
   "source": [
    "unique(objH$metadata$celltype.coarse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objH <- readRDS(\"./cache/CosMxcoarseGrainharmonyObj_markers.RDS\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objH$metadata$celltype.coarse <- as.ordered(factor(objH$metadata$celltype.coarse, c(\"Adipocytes\",\"Fibroblasts\",\"Mural\",\"Endothelial\",\n",
    "                                                                                    \"B cells\",\"Plasma\",\"T cells\",\"ProliferationT\",\"Mast\",\"Myeloid\",\n",
    "                                                                                    \"DC\")))\n",
    "                                                                    \n",
    "cols <- c(\"#ffeca9\",\"gold\",\"sienna1\",\"red3\",\"pink\",\"deeppink\",\"green3\",\"green4\",\"lightblue1\",\"royalblue1\",\"royalblue4\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc9eeeb1-701e-4dff-92c1-5ffde67d5349",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(7, 9)\n",
    "umap_after<-plot_dim_red(\n",
    "    dim_red_embeddings = objH$Humap$embedding,\n",
    "    clusters = objH$Humap$clusters,\n",
    "    metadata = objH$metadata,\n",
    "    cell_id_colname = \"cellID\",\n",
    "    color_by = \"celltype.coarse\",\n",
    "    plot_title = \"UMAP after Harmony\",\n",
    "    dim_red_type = \"UMAP\",\n",
    "    size_points = 0.2,\n",
    "    plot_labels = FALSE\n",
    ")\n",
    "umap_after + theme_powerpoint() + scale_fill_manual(values=cols)+ scale_color_manual(values=cols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c98c36d5-0e73-42e2-831e-1fc9d5b38ce1",
   "metadata": {},
   "outputs": [],
   "source": [
    "dim(objH$counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be97dbaf-3ea0-47a6-883e-8b230c8efb49",
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of cells of each type\n",
    "table(objH$metadata$celltype.coarse) %>% as.data.frame %>% arrange(-Freq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(objH$metadata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objH$metadata$SampleID <- objH$metadata$SampleID.x\n",
    "objH$metadata$FOV <- objH$metadata$FOV.x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9913707f-0a61-421a-b045-187139009db3",
   "metadata": {},
   "source": [
    "### glmer markers for annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "presto.presto <- function (formula, design, response, size_varname, features = NULL, \n",
    "    effects_cov = c(\"\"), ncore = 1, nsim = 100, family = \"poisson\", \n",
    "    min_sigma = 0, verbose = 0L) \n",
    "{\n",
    "    if (is.null(features)) {\n",
    "        features <- rownames(response)\n",
    "    }\n",
    "    if (family %in% c(\"poisson\", \"binomial\", \"nb\")) {\n",
    "        message(\"CAUTION: if using GLMM, make sure your counts are integers!\")\n",
    "    }\n",
    "    design$EXPOSURE <- design[[size_varname]]\n",
    "    fstr <- gsub(size_varname, \"EXPOSURE\", as.character(formula))\n",
    "    formula <- as.formula(sprintf(\"%s~%s\", fstr[[2]], fstr[[3]]), \n",
    "        env = .GlobalEnv)\n",
    "    if (verbose > 0) {\n",
    "        message(\"Set up models\")\n",
    "    }\n",
    "    model_base <- fit_model.presto(formula, design, response[features[[1]], \n",
    "        ], family)\n",
    "    priornames_df <- as.data.frame(VarCorr(model_base))[, 1:3]\n",
    "    if (isGLMM(model_base)) {\n",
    "        priornames_df <- rbind(priornames_df, tibble(grp = \"Residual\", \n",
    "            var1 = NA, var2 = NA))\n",
    "    }\n",
    "    has_offset <- !all(map_lgl(model_base@resp$offset, identical, \n",
    "        0))\n",
    "    betanames_df <- make_betanames_df(model_base, has_offset)\n",
    "    features <- intersect(features, rownames(response))\n",
    "    if (ncore == 1) {\n",
    "        future::plan(sequential)\n",
    "    }\n",
    "    else if (ncore %in% c(0, Inf)) {\n",
    "        ncore <- availableCores()\n",
    "        future::plan(multisession)\n",
    "    }\n",
    "    else {\n",
    "        .ncore <<- ncore\n",
    "        future::plan(future::multisession(workers = .ncore))\n",
    "        rm(.ncore)\n",
    "    }\n",
    "    if (verbose > 0) {\n",
    "        message(\"Learn the models\")\n",
    "    }\n",
    "    lres <- furrr::future_map(features, glmm_uni, formula, design, \n",
    "        response, effects_cov, family, nsim, has_offset, min_sigma)\n",
    "    names(lres) <- features\n",
    "    lres <- lres[which(purrr::map_lgl(as.integer(map_int(lres, \n",
    "        \"status\")), identical, 0L))]\n",
    "    if (verbose > 0) {\n",
    "        message(\"Aggregate the results\")\n",
    "    }\n",
    "    res <- collapse_lres(lres)\n",
    "    if (verbose > 0) {\n",
    "        message(\"Cleap up names\")\n",
    "    }\n",
    "    covmat_names <- tibble(grpvar_orig = rownames(res$covmat)) %>% \n",
    "        left_join(subset(betanames_df, term %in% c(\"(Intercept)\", \n",
    "            \"Fixed\"))) %>% dplyr::mutate(newname = case_when(is.na(grpvar) ~ \n",
    "        grpvar_orig, TRUE ~ as.character(glue::glue(\"{grpvar}.{grp}.{term}\")))) %>% \n",
    "        with(newname)\n",
    "    dimnames(res$covmat) <- list(covmat_names, covmat_names, \n",
    "        colnames(res$beta))\n",
    "    res$betanames_df <- betanames_df\n",
    "    res$priornames_df <- priornames_df\n",
    "    res$meta_data <- design\n",
    "    if (has_offset) {\n",
    "        res$design <- list(EXPOSURE = model_base@resp$offset, \n",
    "            t(model_base@pp$X), model_base@pp$Zt) %>% purrr::reduce(Matrix::rbind2)\n",
    "    }\n",
    "    else {\n",
    "        res$design <- list(t(model_base@pp$X), model_base@pp$Zt) %>% \n",
    "            purrr::reduce(Matrix::rbind2)\n",
    "    }\n",
    "    row.names(res$design) <- res$betanames_df$grp\n",
    "    res$response <- response[names(lres), ]\n",
    "    if (verbose > 0) {\n",
    "        message(\"Compute gene means\")\n",
    "    }\n",
    "    res <- genemeans.presto(res, xpm = 1e+06)\n",
    "    res$has_offset <- has_offset\n",
    "    res$family <- family\n",
    "    res$size_varname <- size_varname\n",
    "    res$nsim <- nsim\n",
    "    res$formula_str <- as.character(formula)\n",
    "    return(res)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#These are internal functions i.e. unexported so we must access with :::\n",
    "collapse_vecs <- presto:::collapse_vecs\n",
    "collapse_mats <- presto:::collapse_mats\n",
    "collapse_lres <- presto:::collapse_lres"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0f27f10-9141-49de-bc14-fb3bf431a7ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "formula_glmer_input<-createFormula(\"y\", \"1 + (1|celltype.coarse)\", c(\"(1|SampleFOV/celltype.coarse)\", \"(1|SampleID/celltype.coarse)\", \"offset(logUMI)\"))\n",
    "\n",
    "objH$markers_celltype_coarse<-find_glmer_markers(\n",
    "    objH$counts,\n",
    "    objH$metadata,\n",
    "    formula_glmer_input, \n",
    "    \"celltype.coarse\", \n",
    "    c(\"SampleID\", \"FOV\", \"celltype.coarse\") \n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d4c0780-3eca-4983-a08b-0381a3341993",
   "metadata": {},
   "outputs": [],
   "source": [
    "objH$markers_celltype_coarse %>% \n",
    "    group_by(contrast) %>% \n",
    "    filter(pvalue < 1e-02) %>% \n",
    "    filter(logFC > 0.5) %>% \n",
    "    top_n(n = 20, wt = logFC) %>% \n",
    "    mutate(rank = rank(-logFC)) %>% \n",
    "    ungroup %>% \n",
    "    # dplyr::mutate(feature = paste0(feature, \"_\", logFC)) %>% \n",
    "    dplyr::select(contrast, feature, rank) %>% \n",
    "    spread(contrast, feature, fill = NA)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90f29693-bd0c-4b20-9c25-7daf1ff9bc49",
   "metadata": {},
   "source": [
    "### heatmap of markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5dec253b-d773-4e50-a550-e566d86f83f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# pick signficant genes to show in the heatmap\n",
    "# arrange the genes such that the gene with max logFC is at the center and values decrease as we move away from center\n",
    "genes_use_df<-data.table(objH$markers_celltype_coarse)[\n",
    "    pvalue < 0.05][\n",
    "    , .SD[order(-zscore)][1, ], by = feature][\n",
    "    , .SD[order(-logFC)][1:20, ], by = contrast][\n",
    "    , .SD[!is.na(feature)]][\n",
    "    , .SD[order(logFC)], by = contrast][\n",
    "    , \":=\" (len = .N, rank = frank(logFC)), by = contrast][\n",
    "    , .SD[c(seq(1, round(.SD$len[1]/3), 1), seq(.SD$len[1], round(.SD$len[1]/3)+1, -1)), ], by = contrast][\n",
    "    , .(feature, logFC, zscore, contrast, pvalue) # the genes in this list will be unique because I assigned one z-score to each gene\n",
    "]\n",
    "\n",
    "message(\"Checking that the genes are unique: \", nrow(genes_use_df) == length(unique(genes_use_df$feature)))\n",
    "message(\"Number of genes in heatmap: \", nrow(genes_use_df))\n",
    "\n",
    "# select the top 3 genes with highest logFC in each cluster and assign each gene to a cell type based on it logFC - the gene belongs to the celltype with higher logFC\n",
    "genes_to_mark_celltypes<-genes_use_df[, .SD[order(-logFC)][1:5, ], by = contrast][\n",
    "    , .(contrast, feature)\n",
    "]\n",
    "genes_to_mark_celltypes<-genes_to_mark_celltypes %>% \n",
    "    # subset(! contrast %in% c(\"Neutrophils\", \"T cells\")) %>% \n",
    "    rbind(genes_use_df[feature %in% c(\"MPO\", \"CYSTM1\", \"EFNA4\"), .(contrast, feature)]) %>% \n",
    "    rbind(genes_use_df[feature %in% c(\"CD3G\", \"CD3E\", \"IL32\"), .(contrast, feature)])\n",
    "\n",
    ".mat<-objH$markers_celltype_coarse %>% \n",
    "    subset(feature %in% genes_use_df$feature) %>% \n",
    "    dplyr::select(cluster, logFC, feature) %>% \n",
    "    tidyr::spread(cluster, logFC) %>% \n",
    "    tibble::column_to_rownames('feature') %>% \n",
    "    as.matrix() %>% t() \n",
    ".mat<-.mat[, genes_use_df$feature]\n",
    "\n",
    "message(\n",
    "    \"Are all genes that are being marked also present in the matrix?: \", \n",
    "    length(setdiff(genes_to_mark_celltypes$feature, colnames(.mat))) == 0\n",
    ")\n",
    "\n",
    "genes_idx<-data.frame(idx = which(colnames(.mat) %in% genes_to_mark_celltypes$feature)) %>% \n",
    "    mutate(feature = colnames(.mat)[idx])\n",
    "\n",
    "genes_to_mark_celltypes<-genes_to_mark_celltypes %>% left_join(genes_idx, by = \"feature\")\n",
    "\n",
    "\n",
    "names(cols)<-char(levels(objH$metadata$celltype.coarse))\n",
    "message(\"Colors being used: \")\n",
    "print(cols)\n",
    "message(\"Check if any cell types above are not getting a color because of a smaller palette\")\n",
    "\n",
    "\n",
    "fig.size(5, 15)\n",
    "ha<-columnAnnotation(\n",
    "    col_ann = anno_mark(\n",
    "        at = genes_to_mark_celltypes$idx,\n",
    "        side = \"bottom\",\n",
    "        labels = genes_to_mark_celltypes$feature)\n",
    "    )\n",
    "\n",
    "ha2<-columnAnnotation(show_legend = FALSE,\n",
    "    celltypes = genes_use_df$contrast,\n",
    "    col = list(celltypes = cols)\n",
    "    )\n",
    "\n",
    "# ha2<-columnAnnotation(\n",
    "#     celltypes = anno_block(\n",
    "#         gp = gpar(fill = col_pal), \n",
    "#         labels = genes_use_df$contrast %>% unique\n",
    "\n",
    "# ))\n",
    "split<-genes_use_df$contrast\n",
    "colnames(.mat)<-NULL\n",
    ".mat %>% \n",
    "    Heatmap(\n",
    "        name = 'logFC',\n",
    "        col = circlize::colorRamp2(c(0, 0.1, 3), c('white', 'white', muted('blue'))),\n",
    "        # column_names_gp = grid::gpar(fontsize = 10),\n",
    "        row_dend_side = \"right\",\n",
    "        cluster_rows = FALSE,\n",
    "        cluster_columns = FALSE,\n",
    "        top_annotation = ha2, \n",
    "        # row_km = 12,\n",
    "        bottom_annotation = ha,\n",
    "        column_split = split\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "char(levels(objH$metadata$celltype.coarse))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c904d0e5-4c2a-4467-8b60-bd7ff6da2501",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(5, 25)\n",
    ".mat<-objH$markers_celltype_coarse %>% \n",
    "    subset(feature %in% genes_use_df$feature) %>% \n",
    "    dplyr::select(cluster, logFC, feature) %>% \n",
    "    tidyr::spread(cluster, logFC) %>% \n",
    "    tibble::column_to_rownames('feature') %>% \n",
    "    as.matrix() %>% t() \n",
    ".mat<-.mat[, genes_use_df$feature]\n",
    ".mat %>% \n",
    "    Heatmap(\n",
    "        name = 'logFC',\n",
    "        col = circlize::colorRamp2(c(0, 0.1, 3), c('white', 'white', muted('blue'))),\n",
    "        column_names_gp = grid::gpar(fontsize = 6),\n",
    "        row_dend_side = \"right\",\n",
    "        cluster_rows = FALSE,\n",
    "        cluster_columns = FALSE,\n",
    "        top_annotation = ha2, \n",
    "        # row_km = 12,\n",
    "        # column_title = \"Heatmap of markers for coarse-grain cell types\", \n",
    "        column_split = split\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89cf6616-9ffd-4c28-984e-01ba81b17e37",
   "metadata": {},
   "outputs": [],
   "source": [
    "names(objH)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea13175f-317a-4025-85f6-187deb43d319",
   "metadata": {},
   "source": [
    "## cache"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9825c4f7-59a3-4c02-96e9-bc484f290ad6",
   "metadata": {},
   "outputs": [],
   "source": [
    "saveRDS(objH, \"./cache/CosMxcoarseGrainharmonyObj_markers.RDS\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31730271-3c79-41e8-8077-0193490dd8d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "objH<-readRDS(\"./cache/CosMxcoarseGrainharmonyObj_markers.RDS\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objH$markers_celltype_coarse"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99a51c3e-2fe2-4af5-b223-87bcd9c628e1",
   "metadata": {},
   "source": [
    "## glmer heatmap with zscores instead of logFC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fd2ac34-4193-4e34-a899-c5e0475c3ea3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# pick signficant genes to show in the heatmap\n",
    "# arrange the genes such that the gene with max logFC is at the center and values decrease as we move away from center\n",
    "genes_use_df<-data.table(objH$markers_celltype_coarse)[\n",
    "    pvalue < 0.05][\n",
    "    , .SD[order(-zscore)][1, ], by = feature][\n",
    "    , .SD[order(-logFC)][1:20, ], by = contrast][\n",
    "    , .SD[!is.na(feature)]][\n",
    "    , .SD[order(logFC)], by = contrast][\n",
    "    , \":=\" (len = .N, rank = frank(logFC)), by = contrast][\n",
    "    , .SD[c(seq(1, round(.SD$len[1]/3), 1), seq(.SD$len[1], round(.SD$len[1]/3)+1, -1)), ], by = contrast][\n",
    "    , .(feature, logFC, zscore, contrast, pvalue) # the genes in this list will be unique because I assigned one z-score to each gene\n",
    "]\n",
    "\n",
    "message(\"Checking that the genes are unique: \", nrow(genes_use_df) == length(unique(genes_use_df$feature)))\n",
    "message(\"Number of genes in heatmap: \", nrow(genes_use_df))\n",
    "\n",
    "# select the top 3 genes with highest logFC in each cluster and assign each gene to a cell type based on it logFC - the gene belongs to the celltype with higher logFC\n",
    "genes_to_mark_celltypes<-genes_use_df[, .SD[order(-logFC)][1:3, ], by = contrast][\n",
    "    , .(contrast, feature)\n",
    "]\n",
    "\n",
    ".mat<-objH$markers_celltype_coarse %>% \n",
    "    subset(feature %in% genes_use_df$feature) %>% \n",
    "    dplyr::select(cluster, zscore, feature) %>% \n",
    "    tidyr::spread(cluster, zscore) %>% \n",
    "    tibble::column_to_rownames('feature') %>% \n",
    "    as.matrix() %>% t() \n",
    ".mat<-.mat[, genes_use_df$feature]\n",
    "\n",
    "message(\n",
    "    \"Are all genes that are being marked also present in the matrix?: \", \n",
    "    length(setdiff(genes_to_mark_celltypes$feature, colnames(.mat))) == 0\n",
    ")\n",
    "\n",
    "genes_idx<-data.frame(idx = which(colnames(.mat) %in% genes_to_mark_celltypes$feature)) %>% \n",
    "    mutate(feature = colnames(.mat)[idx])\n",
    "\n",
    "genes_to_mark_celltypes<-genes_to_mark_celltypes %>% left_join(genes_idx, by = \"feature\")\n",
    "\n",
    "message(\"Colors being used: \")\n",
    "message(\"Check if any cell types above are not getting a color because of a smaller palette\")\n",
    "\n",
    "\n",
    "fig.size(5, 15)\n",
    "ha<-columnAnnotation(\n",
    "    col_ann = anno_mark(\n",
    "        at = genes_to_mark_celltypes$idx,\n",
    "        labels = genes_to_mark_celltypes$feature,\n",
    "        side = \"bottom\"\n",
    "        )\n",
    "    )\n",
    "\n",
    "ha2<-columnAnnotation(\n",
    "    celltypes = genes_use_df$contrast,\n",
    "    col = list(celltypes = cols)\n",
    "    )\n",
    "\n",
    "# ha2<-columnAnnotation(\n",
    "#     celltypes = anno_block(\n",
    "#         gp = gpar(fill = col_pal), \n",
    "#         labels = genes_use_df$contrast %>% unique\n",
    "\n",
    "# ))\n",
    "split<-genes_use_df$contrast\n",
    "colnames(.mat)<-NULL\n",
    ".mat %>% \n",
    "    Heatmap(\n",
    "        name = \"zscore\",\n",
    "        col = circlize::colorRamp2(c(-20, 1, 20), c('white', 'white', muted('blue'))),\n",
    "        # column_names_gp = grid::gpar(fontsize = 10),\n",
    "        row_dend_side = \"right\",\n",
    "        cluster_rows = FALSE,\n",
    "        cluster_columns = FALSE,\n",
    "        top_annotation = ha2, \n",
    "        # row_km = 12,\n",
    "        # column_title = \"Heatmap of markers for coarse-grain cell types\", \n",
    "        bottom_annotation = ha,\n",
    "        column_split = split\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2be0220-813b-4c24-879e-9ab011af4167",
   "metadata": {},
   "source": [
    "# Markers in UMAP space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adca4fd0-514f-4e6b-bceb-e44c63243d1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "unique(objH$metadata$celltype.coarse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f41cfd9-7c1c-4d7b-8da3-1a1eee75ffbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "genes_viz<-genes_to_mark_celltypes$feature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f420b7ed-9654-4c8c-b0ec-508f8154f646",
   "metadata": {},
   "outputs": [],
   "source": [
    "cellids<-objH$metadata %>% subset(!grepl(\"Doublet\", celltype.coarse)) %>% with(cellID)\n",
    "counts_metadata<-objH$logcpx[\n",
    "    genes_viz, cellids] %>% \n",
    "    t() %>% \n",
    "    as.matrix %>% \n",
    "    as.data.frame %>% \n",
    "    rownames_to_column(\"cellID\") %>% \n",
    "    gather(target, geneval, -cellID) %>% \n",
    "    left_join(\n",
    "        genes_to_mark_celltypes %>% dplyr::select(contrast, feature), \n",
    "        by = c(\"target\" = \"feature\")\n",
    "    ) %>%  \n",
    "    mutate(celltypes = paste0(target, \" (\", contrast, \")\")) %>% \n",
    "    identity()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d6f3593-aab7-4c76-94e9-849187ab7cd2",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "head(counts_metadata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "markers_umap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "418b50dd-ea06-4505-9f39-0c8a12347823",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(10, 10)\n",
    "p<-markers_umap(\n",
    "    dim_red_embeddings = objH$Humap$embedding[cellids, ],\n",
    "    metadata = counts_metadata, \n",
    "    cell_id_colname = \"cellID\", \n",
    "    color_by = \"geneval\", \n",
    "    plot_title = \"Markers in UMAP space\",\n",
    "    dim_red_type = \"UMAP\"\n",
    ") "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14286576-414b-41bc-8974-51a5ae5134f5",
   "metadata": {},
   "source": [
    "## order the facets based on cell types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35c6736b-0331-4b12-8c26-bf1bb4e17968",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(12, 15)\n",
    "p + facet_wrap(~ contrast + target, labeller = labeller(celltypes = (label_wrap_gen(10))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "genes <- as.data.frame(rownames(objH$counts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.matrix.max.rows=1000, repr.matrix.max.cols=200)\n",
    "genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf06cacc-925a-46e5-b6b5-53dec1c65189",
   "metadata": {},
   "outputs": [],
   "source": [
    "genes_viz<-c(\"IL15\",\"IL2RA\",\"IL2RB\",\"IL2RG\",\"IL15RA\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "objH$metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dimlist(objH)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objH$logcpx[\n",
    "    genes_viz, colnames(objH$logcpx) %in% objH$metadata$cellID\n",
    "    ] %>% \n",
    "    t() %>% \n",
    "    as.matrix %>% \n",
    "    as.data.frame %>% \n",
    "    rownames_to_column(\"cellID\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b815e4d8-5d6f-47a0-8e8b-951b0594718a",
   "metadata": {},
   "outputs": [],
   "source": [
    "counts_metadata<-objH$logcpx[\n",
    "    genes_viz, colnames(objH$logcpx) %in% objH$metadata$cellID\n",
    "    ] %>% \n",
    "    t() %>% \n",
    "    as.matrix %>% \n",
    "    as.data.frame %>% \n",
    "    rownames_to_column(\"cellID\") %>% \n",
    "    gather(target, geneval, -cellID) %>% \n",
    "    identity() %>% arrange(geneval) %>% left_join(objH$metadata %>% dplyr::select(cellID, celltype.coarse), by = \"cellID\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts_metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(10, 20)\n",
    "ggplot(counts_metadata, aes(x=celltype.coarse, y=geneval)) + \n",
    "  geom_violin() + facet_grid(~target) + theme(axis.text.x=element_text(angle=90))  +stat_summary(fun.y=mean, geom=\"point\", size=2, color=\"red\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "markers_umap <- function (dim_red_embeddings, metadata, cell_id_colname, color_by, \n",
    "    label_by, plot_title, dim_red_type, legend_posn = \"right\", \n",
    "    shape_points = \".\", size_points = 0.1) \n",
    "{\n",
    "    plt_df <- dim_red_embeddings[, 1:2] %>% as.data.frame %>% \n",
    "        purrr::set_names(\"V1\", \"V2\") %>% rownames_to_column(cell_id_colname) %>% \n",
    "        left_join(metadata %>% dplyr::select(-starts_with(\"V\")), \n",
    "            by = cell_id_colname, multiple = \"all\") %>% sample_frac(1) %>% arrange(geneval)\n",
    "    range_col <- plt_df %>% with(geneval) %>% range\n",
    "    p <- plt_df %>% ggplot() + geom_point(aes(V1, V2, color = geneval), \n",
    "        shape = shape_points) + ggtitle(plot_title) + scale_color_gradient2_tableau(limits = c(-range_col[2], \n",
    "        range_col[2])) + xlab(paste0(dim_red_type, \"_1\")) + ylab(paste0(dim_red_type, \n",
    "        \"_2\")) + theme(plot.title = element_textbox_simple()) + \n",
    "        theme(legend.position = legend_posn) + facet_wrap(~target) + \n",
    "        NULL\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f0260ca-270c-4128-a07c-f584b864bea6",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(10, 10)\n",
    "p_adip<-markers_umap(\n",
    "    dim_red_embeddings = objH$Humap$embedding,\n",
    "    metadata = counts_metadata, \n",
    "    cell_id_colname = \"cellID\", \n",
    "    color_by = \"geneval\", \n",
    "    plot_title = \"Markers in UMAP space\",\n",
    "    dim_red_type = \"UMAP\"\n",
    ") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c37e374-45b3-4874-b834-848b11a18aed",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(7, 15)\n",
    "p_adip"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c125c84-415a-4370-8893-8d7f6e8078d1",
   "metadata": {},
   "source": [
    "# Coarse cell types in space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cellgeoms<-readRDS(\"/datastore/lucy/CosMx/cache/raw/Glasgow_cellgeoms_QC.RDS\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cellgeoms<-cellgeoms %>%\n",
    "    left_join(objH$metadata %>% dplyr::select(cellID, celltype.coarse), by = \"cellID\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cellgeoms <- cellgeoms %>% \n",
    "    unite(\"SampleFOV\", c(\"SampleID\", \"FOV\"), remove = FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cellgeoms <- merge(cellgeoms, sampleFOVcondition, all=FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cellgeoms$SampleFOV <- paste0(cellgeoms$Condition,\"_\",cellgeoms$SampleFOV)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37a2dfe5-6886-4c84-b095-2a0daed1e6cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "unique(cellgeoms$celltype.coarse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(cellgeoms %>% subset(is.na(celltype.coarse)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "names(cols) <- levels(cellgeoms$cell.coarse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da6191f7-7ffd-4c79-8fbd-38076ac775ab",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "library(ggdark)\n",
    "options(repr.plot.width=6, repr.plot.height=6)\n",
    "cellgeoms %>% \n",
    "    with(SampleFOV) %>%  \n",
    "    unique %>% \n",
    "    # head(1) %>% \n",
    "    purrr::map(~ cellgeoms %>% \n",
    "        # subset(!grepl(\"Doublet*|Low_Quality\", celltype.coarse)) %>% \n",
    "        subset(SampleFOV %in% .x) %>% \n",
    "        ggplot() + \n",
    "            geom_sf(aes(geometry = shape, fill = celltype.coarse), na.rm=TRUE,linewidth = 0.1) + \n",
    "            scale_fill_manual(values = cols, na.value = \"black\", na.translate = TRUE) + \n",
    "            ggtitle(.x) + dark_theme_classic() \n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d4f6891-9d02-455a-b294-f41a788871dd",
   "metadata": {},
   "source": [
    "# Cell fractions across samples and FOVs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampleFOVcondition <- read.table(\"/datastore/lucy/CosMx/cache/SampleFOV-condition.csv\", header=TRUE, sep=\",\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objH$metadata <- objH$metadata %>% inner_join(sampleFOVcondition, by = \"SampleFOV\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(8, 4)\n",
    ".vals<-unique(objH$metadata$celltype.coarse)\n",
    ".vals<-.vals[!is.na(.vals)]\n",
    "\n",
    "objH$metadata %>% \n",
    "    with(table(celltype.coarse, Condition)) %>% \n",
    "    prop.table(margin = 2) %>% \n",
    "    as.data.frame %>% \n",
    "    ggplot() + \n",
    "        geom_col(aes(x = Condition, y = Freq, fill = celltype.coarse)) + \n",
    "        scale_fill_manual(values = cols) +\n",
    "        guides(x =  guide_axis(angle = 90)) + theme_powerpoint() + theme(axis.text=element_text(size=20), axis.title=element_text(size=30))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c153e32f-a262-4e61-813f-7d8b10053cd8",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(5, 5)\n",
    "\n",
    "cellgeoms %>% \n",
    "    unite(\"SampleFOV\", c(\"SampleID\", \"FOV\"), remove = FALSE) %>% \n",
    "    with(table(SampleID, SampleFOV, celltype.coarse)) %>% \n",
    "    prop.table(margin = 1) %>% \n",
    "    as.data.frame %>% \n",
    "    ggplot() + \n",
    "        geom_col(aes(SampleID, Freq, fill = celltype.coarse)) + \n",
    "        scale_fill_manual(values = cols) + \n",
    "        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3aae9de7-2dc0-4969-9a22-d5f63c282a3c",
   "metadata": {},
   "source": [
    "## cell fractions across FOVs in each of the three slides"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3abe894-289d-4a94-ae69-aa76402b0944",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(10, 10)\n",
    "cellgeoms %>% \n",
    "    unite(\"SampleFOV\", c(\"SampleID\", \"FOV\"), remove = FALSE) %>% \n",
    "    subset(SampleID %in% \"Run5654_399G0_32G77\") %>% \n",
    "    with(table(SampleID, SampleFOV, celltype.coarse)) %>% \n",
    "    prop.table(margin = 2) %>% \n",
    "    as.data.frame %>% \n",
    "    ggplot() + \n",
    "        geom_col(aes(SampleFOV, Freq, fill = celltype.coarse)) + \n",
    "        scale_fill_manual(values = cols) + \n",
    "        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + \n",
    "        ggtitle(\"Fraction of cells in each FOV\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb72320b-fb9b-4b1f-b4aa-a37781f4425d",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(10, 10)\n",
    "cellgeoms %>% \n",
    "    unite(\"SampleFOV\", c(\"SampleID\", \"FOV\"), remove = FALSE) %>% \n",
    "    subset(SampleID %in% \"Run5666_G2697_G1532\") %>% \n",
    "    with(table(SampleID, SampleFOV, celltype.coarse)) %>% \n",
    "    prop.table(margin = 2) %>% \n",
    "    as.data.frame %>% \n",
    "    ggplot() + \n",
    "        geom_col(aes(SampleFOV, Freq, fill = celltype.coarse)) + \n",
    "        scale_fill_manual(values = cols) + \n",
    "        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + \n",
    "        ggtitle(\"Fraction of cells in each FOV\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1497f659-bfac-41e9-b316-ccd06375955c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.size(10, 10)\n",
    "cellgeoms %>% \n",
    "    unite(\"SampleFOV\", c(\"SampleID\", \"FOV\"), remove = FALSE) %>% \n",
    "    subset(SampleID %in% \"Run5666_G7G72_G22G5\") %>% \n",
    "    with(table(SampleID, SampleFOV, celltype.coarse)) %>% \n",
    "    prop.table(margin = 2) %>% \n",
    "    as.data.frame %>% \n",
    "    ggplot() + \n",
    "        geom_col(aes(SampleFOV, Freq, fill = celltype.coarse)) + \n",
    "        scale_fill_manual(values = cols) + \n",
    "        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + \n",
    "        ggtitle(\"Fraction of cells in each FOV\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add in finetyping results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "integrated.tcell <- readRDS(\"./cache/integrated.tcell.rds\")\n",
    "integrated.fibro <- readRDS(\"./cache/integrated.fibro.RDS\")\n",
    "integrated.myeloid <- readRDS(\"./cache/integrated.myeloid.RDS\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "integrated.fibro$metadata$type_harmony <- integrated.fibro$metadata$type_harmony_clean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "integrated.tcell$metadata$type_harmony <- paste0(\"Tcell_\",integrated.tcell$metadata$type_harmony)\n",
    "integrated.myeloid$metadata$type_harmony <- paste0(\"Myeloid_\",integrated.myeloid$metadata$type_harmony)\n",
    "integrated.fibro$metadata$type_harmony <- paste0(\"Fibro_\",integrated.fibro$metadata$type_harmony)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(integrated.tcell$metadata)\n",
    "head(integrated.fibro$metadata)\n",
    "head(integrated.myeloid$metadata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(objH$metadata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fine.meta <- objH$metadata %>% left_join(\n",
    "    rbind(integrated.tcell$metadata %>% dplyr::select(cellID, type_harmony),\n",
    "          integrated.myeloid$metadata %>% dplyr::select(cellID, type_harmony),\n",
    "          integrated.fibro$metadata %>% dplyr::select(cellID, type_harmony)\n",
    "    ), \n",
    "    by = \"cellID\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(fine.meta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objH$metadata <- fine.meta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Replace NA values for celltypes not finetyped\n",
    "objH$metadata$type_harmony <- ifelse(is.na(objH$metadata$type_harmony), objH$metadata$celltype.coarse, objH$metadata$type_harmony)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "saveRDS(objH, \"./cache/CosMxcoarseGrainharmonyObj_markers.RDS\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objH <- readRDS(\"./cache/CosMxcoarseGrainharmonyObj_markers.RDS\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objH$metadata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fine cell types in space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(sf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cellgeoms<-readRDS(\"/datastore/lucy/CosMx/cache/raw/Glasgow_cellgeoms_QC.RDS\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(cellgeoms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cellgeoms<-cellgeoms %>%\n",
    "    left_join(objH$metadata %>% dplyr::select(cellID, type_harmony, celltype.coarse), by = \"cellID\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampleFOVcondition <- read.table(\"/datastore/lucy/CosMx/cache/SampleFOV-condition.csv\", header=TRUE, sep=\",\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cellgeoms <- cellgeoms %>% \n",
    "    unite(\"SampleFOV\", c(\"SampleID\", \"FOV\"), remove = FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cellgeoms <- merge(cellgeoms, sampleFOVcondition, all=FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cellgeoms$SampleFOV <- paste0(cellgeoms$Condition,\"_\",cellgeoms$SampleFOV)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique(cellgeoms$type_harmony_clean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusters.of.interest <- c(\"Myeloid_CD63+ DC3\",\"Myeloid_CLEC10A+ STM\",\"Myeloid_CLEC9A+ DC1\",\"Myeloid_DC2\",\"Myeloid_CCR7+ migDC\",\n",
    "                  \"Tcell_CCL5+CXCR6+MAF+ TPH\",\"Tcell_MAF+ TPH\",\"Tcell_RORA+CD4+ TEM\",\"Tcell_HSP+CD4+ TEM\",\n",
    "                 \"Tcell_ANXA1+CD4+ TEM\",\"Tcell_CCR7+LTB+CD4+ TCM\",\"Tcell_CCR7+SELL+CD4+ Naive\",\"Tcell_FOXP3+ Treg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cellgeoms$type_dct <- ifelse(cellgeoms$type_harmony==clusters.of.interest, cellgeoms$type_harmony, NA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num.clust <- length(unique(cellgeoms$type_dct))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cols <- tableau_color_pal('Tableau 20')(num.clust)\n",
    "cols <- rainbow(n=num.clust)\n",
    "names(cols) <- unique(cellgeoms$type_dct)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "library(ggdark)\n",
    "options(repr.plot.width=10, repr.plot.height=10)\n",
    "cellgeoms %>% \n",
    "    with(SampleFOV) %>%  \n",
    "    unique %>% \n",
    "    # head(1) %>% \n",
    "    purrr::map(~ cellgeoms %>% \n",
    "        subset(SampleFOV %in% .x) %>% \n",
    "        ggplot() + \n",
    "            geom_sf(aes(geometry = shape, fill = type_dct), na.rm=TRUE,linewidth = 0.1) + \n",
    "            scale_fill_manual(values = cols, na.value = \"black\", na.translate = TRUE) + \n",
    "            ggtitle(.x) + dark_theme_classic() \n",
    "    )"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.1.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}